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Series Foreword 


Dear reader, 

The series Simula SpringerBriefs on Computing was established in 2016, with 
the aim of publishing compact introductions and state-of-the-art overviews of se- 
lect fields in computing. Research is increasingly interdisciplinary, and students and 
experienced researchers both often face the need to learn the foundations, tools, 
and methods of a new field. This process can be demanding, and typically involves 
extensive reading of multidisciplinary publications with different notations, termi- 
nologies and styles of presentation. The briefs in this series are meant to ease the 
process by explaining important concepts and theories in a specific interdisciplinary 
field without assuming extensive disciplinary knowledge and by outlining open re- 
search challenges and posing critical questions in the field. 

Simula has a major research program in computational physiology that includes 
a long and close collaboration with the University of California (UC) San Diego. To 
reflect this research focus, we established in 2020 a new subseries entitled Simula 
Springer Briefs on Computing - Reports on Computational Physiology. The subseries 
includes both introductory and advanced texts on select fields of computational 
physiology, designed to advance interdisciplinary scientific literacy and promote 
effective communication and collaboration in the field. This subseries is also the 
outlet for collections of reports from the annual Summer School in Computational 
Physiology, organized by Simula, University of Oslo, and UC San Diego. The school 
starts in June each year with students spending two weeks in Oslo learning the 
principles underlying mathematical models commonly used in studying the heart 
and the brain. During their stay in Oslo, students are assigned a research project to 
work on over the summer. In August, they travel to San Diego for another week of 
training and project work, and a final presentation of their findings. Every year, we 
have been impressed by the students' creativity and we often see results that could 
lead to a scientific publication. Starting with the 2021 edition of the summer school, 
we have taken the course one step further by having each team conclude their project 
with a scientific report that can pass rigorous peer review as a publication in this 
new series. 


vi Series Foreword 


All items in the main series and the subseries are published within the 
SpringerOpen framework, as this will allow authors to use the series to publish 
an initial version of their manuscript that could subsequently evolve into a full-scale 
book on a broader theme. Since the briefs are freely available online, the authors do 
not receive any direct income from the sales; however, remuneration is provided for 
every completed manuscript. Briefs are written on the basis of an invitation from a 
member of the editorial board. 

Suggestions for possible topics are most welcome, and interested authors are 
encouraged to contact a member of the editorial board. 


March 2023 Dr. Joakim Sundnes 
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Preface 


Since 2014, we have organized an annual summer school in computational physiol- 
ogy. The school starts in June each year and the graduate students spend two weeks 
in Oslo learning the principles underlying mathematical models commonly used in 
studying the heart and the brain. At the end of their stay in Oslo, the students are 
assigned a research project to work on over the summer. In August the students travel 
to the University of California, San Diego to present their findings. Each year, we 
have been duly impressed by the students’ progress and we have often seen that the 
results contain the rudiments of a scientific paper. 

Starting in the 2021 edition of the summer school, we have taken the course one 
step further and aim to conclude every project with a scientific report that passes 
rigorous peer review as a publication in this new series called Simula SpringerBriefs 
on Computing — reports on computational physiology. 

One advantage of this course adjustment is that we have the opportunity to intro- 
duce students to scientific writing. To ensure the students get the best introduction 
in the shortest amount of time, we have commissioned a professional introduction to 
science writing by Nature. The students participate in a Nature Masterclasses work- 
shop in order to strengthen skills in high quality scientific writing and publishing. 
The workshop is tailored to publications in the field of computational physiology 
and allows students to gather real-time feedback on their reports. 

We would like to emphasise that the contributions in this series are brief re- 
ports based on the intensive research projects assigned during the summer school. 
Each report addresses a specific problem of importance in physiology and presents 
a succinct summary of the findings (8-15 pages). We do not require that results 
represent new scientific results; rather, they can reproduce or supplement earlier 
computational studies or experimental findings. The physiological question under 
consideration should be clearly formulated, the mathematical models should be de- 
fined in a manner readable by others at the same level of expertise, and the software 
used should, if possible, be made publicly available. All reports in this series are 
subjected to peer-review by the other students and supervisors in the program. 


vii 
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Chapter 1 updates 
Studying the Role of Astrocytic Membrane 
Properties on Microscopic Fluid Flow in Brain 
Tissue 


Nigar Abbasova, Adyant Balaji, Elena Bernardelli, Ada J Ellingsrud, and Marte J 
Setra 


Abstract The role of astrocyte networks in brain volume homeostasis and waste 
clearance has not received enough attention from the neuroscience community. 
However, recent research efforts indicate that glial cells are crucial for fluid flow 
through brain tissue, contributing to clearance and maintenance of brain volume. 
We examine the role of various glial cotransporters in the spatial and temporal 
changes of the intra- and extracellular volume fractions and fluid dynamics via 
computational modelling. The model is incorporated within the Kirchhoff-Nernst- 
Planck electrodiffusive framework and takes into account ionic electrodiffusion and 
fluid dynamics. Our research shows that all model configurations demonstrate similar 
fluid fluxes, except those involving HCO; dynamics. The model configuration that 
included the NBC cotransporter was observed to have the greatest intracellular total 
volume-weighted fluid velocity of 16 um/s. 
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2 Microscopic Fluid Flow in the Brain 


1.1 Introduction 


Although astrocytes comprise a large part of the glial cell population in the mam- 
malian brain [1], the role of astrocyte networks has not received enough attention 
from the neuroscience community. However, recent research within astrocytic net- 
works and brain clearance pathways shows that astrocytes are likely to be crucial 
for fluid flow through brain tissue, which contributes primarily to brain volume 
homeostasis and waste clearance [2]. One of the key questions that has yet to be 
understood relates to the role of the astrocytic membrane and its effect on the flow 
of transmembrane and compartmental fluid. Another closely related phenomenon is 
the shrinkage of the extracellular space between neurons and surrounding astrocytes 
during neuronal stimulation [3]. 

Ostby et al. 2009 [3] studied the extracellular shrinkage induced by neural activity 
by investigating the astrocytic membrane mechanisms. They report that Na*/K*/CI ^ 
(NKCC1) and Na*/HCO; (NBC) cotransporters appear to be crucial factors in 
achieving the measurable extent of the shrinkage of the extracellular space (ECS). 
Jin et al. [4] studied the dynamics of extracellular K * ions during neural stimulation 
and the effect of the membrane water channel aquaporin-4 (AQP4) in astrocytes. 
However, both authors ([3, 4]) excluded spatial effects from consideration, thus 
preventing them from looking into intracompartmental gradients. Afterwards, Sætra 
et al. 2023 [2] introduced a computational framework to characterise the spatial and 
temporal dynamics of the astrocytic network triggered by neural activity. However, 
the authors were able to estimate the electrochemomechanical response as well as the 
water movement within the compartments of the astrocytic network without studying 
the effects of different membrane mechanisms. 

In this work, we implement other membrane mechanisms in the framework of 
Sætra et al. 2023 [2] and evaluate how they affect microscopic fluid flow in brain 
tissue. We do so by analysing the fluid velocities in the astrocyte network and the ex- 
tracellular space surrounding it. We also investigate the effect of the new membrane 
mechanisms on the spatial and temporal dynamics of the intra- and extracellular vol- 
ume fractions. Our base model for membrane mechanisms includes Na*, K* and CI ^ 
leak channels with a Na*/K* pump. We explore other model configurations by ex- 
panding on the base model by adding different combinations of the K*-coupled C17 
transporter (KCC1), the Na*-coupled HCO; transporter (NBC) and the Na*/K*/CI ^ 
(NKCC1) cotransporter. The different model configurations are adopted from Østby 
et al. 2009 [3]. Our research revealed that all models demonstrate similar ionic fluxes, 
apart from those involving HCO; dynamics. Additionally, these model configura- 
tions had the highest total fluid velocities, both in the intra- and extracellular spaces. 
The highest intracellular and extracellular fluid velocity observed was 16 Lum/min in 
the model configuration that involved the NBC cotransporter. 
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1.2 Methods 
1.2.1 Modelling Fluid Flow Through Astrocyte Networks 


Our aim is to predict the temporal and spatial evolution of the volume fraction o;., 
ion concentrations [Na*],, [K*],, [C17]; and [HCOj ];, electric potential ¢,, and 
hydrostatic pressure p; in the intracellular space (ICS, r = 7) and the extracellular 
space (ECS, r = e). To do so, we build on a model developed by Sætra et al. [2], 
which models the system as a 1D domain of length 300 um representing the tissue 
between two blood vessels. See Figure 1 for a visualisation of the system. 

To model compartmental fluid flow, we implement the M3 model scenario from 
Sætra et al. 2023 [2]. The intracellular fluid flow is driven by hydrostatic and osmotic 
forces, while the extracellular fluid flow is driven by hydrostatic and electro-osmotic 
forces. The fundamental basis of the model is the electrodiffusive Kirchhoff-Nernst- 
Planck framework, and the model dynamics is described using coupled partial dif- 
ferential equations. See Sætra et al. 2023 [2] Section 4.2 for more details. 


;Na ;Na 
Jinput Jdecay 


‘Na Cl i| ; 
Jieak Jleak Jkir Jpump Wm 


ai, [k]i, ĝi, pi j — 


50um 


Fig. 1: Representation of the model devised by Sætra et al. (reproduced with per- 
mission from Sætra et al. 2023) [2]. A: Brain tissue between two blood vessels with 
astrocytes (purple), neurons (grey), and ECS with neuronal activity in the middle. 
B: System represented as 1D domain, including ICS (astrocytes) and ECS. The neu- 
ronal activity is represented by the input currents of Na* and K* ( din and FN in 
the input zone and the decay currents ( dc and dad throughout the domain. The 
transmembrane currents are modelled as an inward rectifying K* current (jxir), Na* 
and Cl” leak currents ( rd and dtu as well as a Na*/K* pump current (pump). 
Intracellular and extracellular currents ( j and JS are driven by electrodiffusion and 
advection. The compartmental fluid flow in the intracellular and extracellular space 
is denoted by uj; and ue, respectively, and the transmembrane fluid flow is denoted 


by Wm. 
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1.2.2 Modelling Neuronal Activity 


Sætra et al. 2023 model neuronal activity as a stimulus in the form of an external 
input flux of K* ions into the ECS and the simultaneous removal of Na* ions from 
active neurons within the input zone. We implement the same method with the input 
zone defined in the interval [L1, L2] with L4 = 1.35: 107^ m and Lo = 1.65- 107^ m, 
during time interval [10 s, 20 s]. We model the input fluxes by following the same 
framework: 


K 2 M 5; 
Jinput = ^ Jinput = Jin- (1.1) 
Here, Jinput is the input flux for potassium, dj is the input flux for sodium and jin 


(mol/ (m?s)) is a constant. 
Similarly, the decay of both of these ions from the extracellular space due to 
neuronal pumps and cotransporters is represented by the decay fluxes, 


ds = ee = —kdec([K*]e — [K" ]ejait); (1.2) 


where the decay current J$ is defined across the whole domain: j$ : Q x 
ecay decay 


(0, T] — R (mol/(n?s)). kaec here denotes the decay factor for the concentration 
of extracellular potassium [K*]., and [K*]. init is the initial extracellular potassium 
concentration. 


1.2.3 Membrane Mechanisms 


In Table 1.1, we present a summary of membrane mechanisms, derived from the 
study by Østby et al. 2009 [3]. We modify the transmembrane flux density of the 
ions j% for each ion species k as described in each model configuration. 


Table 1.1: Model configurations implemented in the study. 


Model configuration Membrane mechanisms 

MCI! (base model) Na*, K+, Cl” leak channels + Na*/K* pump 
MC2 MCI + cotransporter KCC1 

MC3 MCI + cotransporter NBC 

MC4 MCI + cotransporter NKCCI 

MC5 MCI! + cotransporters NBC and NKCC1 


For all the parameters in the following sections, we refer to Table 1.2 in Section 
1.2.4. 
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1.2.3.1 Model Configuration 1 


Model configuration 1 (MC1) models the transport of water by osmosis across the 
membrane, as well as hydrostatic pressure gradients. We also adopt the mechanisms 
of the ionic transport, such as a Na* leak channel, a C17 leak channel, a K* leak 
channel, and a Na*/K* pump. The membrane flux densities (mol/ (m?s)) are given 
by: 


in = a (m — Ewa) + 3 jpump> (1.3) 
ZNa 
du = EX (Om z Ek) ~ 2 jpump: (1.4) 
ZK 
Cl &§Cl 
= =~ (¢m - Eo). 1.5 
Jm Fea (dm ci) (1.5) 
The reversal potentials are given by the Nernst equation: 
RT HS 
k 7 5 ln , (1.6) 
Fz | [k]i 


where R is the gas constant, T is the temperature at which we compute the reversal 
potential and the pump flux density jpump is given by: 


[Na*];^ \ [K*]. | 
[Na*]:? + PS [K*].+ Px, ^ 


Jpump = Pro| (1.7) 


where Ppump is the maximum pump rate, Pya; is the [N a* |; threshold and Px, is the 
[K*]. threshold. 


1.2.3.2 Model Configuration 2 


Model configuration 2 (MC2) is a derivation of MCI, with an additional cotransporter 
KCCI responsible for the transport of K* and CI^ both inward and outward across 
the neuronal membrane. We implement the cotransporter by modifying equations 
(1.4) and (1.5) by adding a term describing the flux of K* and Cl” ions through the 
cotransporter given as 


(1.8) 


RT K+*]e [C17 
dicor oe » w(! Jel 3! 


1.2.3.3 Model Configuration 3 
Model configuration 3 (MC3) is a derivation of MCI, with an additional cotransporter 


NBC responsible for the transport of Na* and HCO; inward across the neuronal 
membrane. 
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We implement the cotransporter by modifying equation (1.3) by adding the term 
describing the flux of Na*: 


SNBC 
( 


JNBC = ———— 
FzNBC 


$m — Esc), (1.9) 


and considering a new equation describing the ion flux of HCO}, that is 


HCO; 
Jm ^ =2Jyngc (1.10) 


Here, gngc is the conductance per unit area for the NBC cotransporter. The reversal 
potential of NBC is 


(1.11) 


RT [Na*].[HCO; H 
EwNBc = ; 


zwNBcF V [Na*];[HCO;3 ]? 
where zwpc is the effective valence of the NBC cotransporter complex, here taken 
to be -1, setting zngc = —(n — 1) = -1 where n is the stoichiometry, and adopting 
n=2., 


1.2.3.4 Model Configuration 4 


Model configuration 4 (MC4) is a derivation of MCI, with an additional cotrans- 
porter NKCCI responsible for the transport of Na*, K* and Cl”. We implement the 
cotransporter by modifying equations (1.3) and (1.4) by adding a term that describes 
the flux of Na* and K* ions, given as 


. &NKkcci RT 


JNKCCI = In | 


ae (1.12) 


We also modify equation (1.5) by adding (1.12) multiplied by a factor of two to 
describe the flux of two Cl ions. Here, gyxcc, is the conductance per unit area for 
the NKCC1 cotransporter. 


1.2.3.5 Model Configuration 5 


Model configuration 5 (MC5) is a derivation of MCI with two additional cotrans- 
porters, NBC and NKCC1, responsible for transport Na*, K*, Cl” and HCO; inward 
across the neuronal membrane. We implement MC5 by modifying equations (1.3), 
(1.4) and (1.5) by adding the two terms describing the flux of Na*, K* and C17, 
given as (1.12) and considering new equations describing the ions flux of HCO; as 
in equations (1.10). 
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1.2.4 Model Implementation 


1.2.4.1 Boundary Conditions 


The boundary conditions incorporated into the model have been adopted from Sætra 
et al. 2023 [2]. This includes sealed-end boundary conditions to ensure that no ions 
or fluids are entering of leaving the system on the boundary I: 


a,j* «np 20 on T, (1.13) 


Q;u, -np =0 on [, (1.14) 


where nr is the outward pointing normal vector. 
We constraint the electrical potential by requiring that 


fo dx — 0, (1.15) 
Q 


where ġo is the extracellular electric potential. To enforce this zero-average con- 
straint, we introduce an additional unknown Lagrange multiplier ce, as described 
by Sætra et al. 2023 in Section 4.7. Following the same set up as in [2] for the 
extracellular hydrostatic pressure pe, we set 


Po=9 on Tight. (1.16) 


1.2.4.2 Initial Conditions 


We first set a collection of pre-established initial values from empirical measurements 
of ion concentrations as shown in Table 1.2 and Table 1.3. We then calibrated the 
models by running simulations for 10^ s, with N = 400 and Ar = 1077, setting 
the transmembrane water permeability 7 to zero. Values of the intracellular and 
extracellular ion concentrations from the final timestep were then used as initial 
conditions throughout the paper. When running the simulations with stimulus, we 
set the water permeability to a nonzero value from Table 1.2, and change the time 
step from At = 107? to At = 107? for a total simulation time of 60 s. 


1.2.4.3 Model Parameters 


We implement the model parameters shown in Table 1.2. 
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Table 1.2: Model parameters. 
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Definition 


Length of domain 

Faraday's constant 

Gas constant 

Temperature 

Na* diffusion constant 

K* diffusion constant 

Cl” diffusion constant 

HCO, diffusion constant 

Na* valence 

K* valence 

Cl” valence 

HCO, valence 

NBC cotransporter valence 
Membrane conductance for Nat 
Membrane conductance for K* 
Membrane conductance for C17 
Membrane conductance for gnxcc1 
Membrane conductance for gxcci 
Membrane conductance for gngc 
Maximum pump rate 

[Na*]; threshold for Na*/K* pump 
[K*]. threshold for Na*/K* pump 
Constant input flux density 

Decay factor for [K* ]. 


Value Ref. 
3.0-107^ m [2] 
96485.3 C/mol 

8.314 J/(molK) 

300K 

1.33-107? m?/s [2] 
1.96-10-? m?/s [2] 
2.03107? m?/s [2] 
1.09-107? m?/s [3] 
1 

1 

-1 

-1 

Si 

1 S/m? [2] 
16.96 S/m? [5] 
1 S/m? [2] 
2-107? S/m? [3] 
7-107! S/m? [3] 
8-107! S/m? [3] 
1.12-1076 mol/(m?s) [2] 
10 mol/ n? [2] 
1.5 mol/ n? [21 
9.05-1077 mol/(m?s) 

2.9.1075 m/s [2] 


The implemented code is openly available at 


https://github.com/hittheant/Simula Summer Project 1. 


1.3 Results 


1.3.1 Calibration 


We calibrate the system as described in the Methods section. Figure 1.2 illustrates 
the temporal and spatial dynamics of the system after calibrating the model and 
running the simulation for 30 s without stimulus for base model MCI. 
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Fig. 1.2: Post-calibration dynamics of the system with model configuration 1 (MC1) 
without stimulus. For this simulation, the water permeability was set to zero. Left two 
columns: Spatial dynamics post-calibration, measured at the end of the simulation 
(t 2 30 s). Right two columns: Time evolution measured at L — 1.50-1074 m. 


We expect the ionic concentrations to remain constant in time and space when 
no stimulus is applied, which is what we observe for MCI (Figure 1.2 C-D). We 
anticipate that our system will demonstrate the same behaviour in space and time 
for other physical quantities, such as variations in the transmembrane hydrostatic 
pressure difference (Figure 1.2 G), membrane potential (Figure 1.2 H) and com- 
partmental volume fractions (Figure 1.2 E-F). The expected result is the same for 
all model setups; however, we only demonstrate the results for MC1 since the other 
model configurations appear to behave the same after calibration. 

The post-calibrated ion concentrations are summarised in Table 1.3 for each model 
configuration. We observe that the ionic concentrations of the different species do not 
vary drastically between the different models. The largest difference appears between 
MC2 and MCS, where we see that MCS has the lowest intracellular Na* concentration 
of 14.873 mM, and MC2 has the highest intracellular Na* concentration of 16.017 
mM. We observe a similar pattern between all ionic species in these two model 
configurations, where MC2 has the highest intracellular concentration of all ions 
and MC5 has the lowest. We also note that MC3 and MC4 have almost identical 
baseline concentrations for all ions, with the exception of small differences of 0.001 
mM between the intracellular K* concentration and the extracellular Na* and K* 
concentrations. 
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Table 1.3: Ionic concentrations after calibration for different model scenarios. Em- 
pirical concentrations are taken from [5] and [6] and were used as initial values for 
the model pre-calibration. 


[Na* ]i [K*]i [CI]; [HCOj] INa'.  [K*]. [Ch]. [HCO}]. 
(mM) (mM) (mM) (mM) (mM) (mM) (mM) (mM) 
Empirical) 15.189 99.959 5.145 11.300 144.662 3.082 133.71 15.000 


MCI 15.510 99.190 5.064 146.580 3.019 133.872 
MC2 16.017 101.812 8.178 140.498 2.844 127.644 
MC3 15.510 99.189 5.064 8.416 146.581 3.020 133.872 14.069 
MC4 15.510 99.190 5.064 146.580 3.019 133.872 


MC5 14.873 98.520 3.691 8.526 149.173 3.040 136.618 13.848 


The membrane potentials are also very similar for all models with $4, = —85 mV, 
with the exception of MC2, where øm = —86 mV (Table 1.4). We also observe the 
same deviation for MC2 and MCS as before, with different jpump values compared 
to the rest of the configurations. Specifically, for MC2, the jpump value is 0.491 
mol/m?s, while MC5 has the lowest Jpump Value of 0.483 mol/m?s. Generally, the 
membrane potential øm, pump flux jpump, and reversal potentials for each ionic 
species appear to remain within a similar range (Table 1.4). We also see on Table 
1.4 that the reversal potential of Cl” varies widely across the different models, with 
the highest and lowest Ec; values 22mV apart. 


Table 1.4: Physical quantities after calibration for different model scenarios. 


om Jpump Ena Ex Eq 

(mV) |(mol/m?s) (mV) (mV)|(mV) 
MCI1| -85 0.493 -58 | 90 | 85 
MC2| -86 0.491 -56 | 92 | 71 
MC3| -85 0.493 -58 | 90 | 85 
MC4| -85 0.493 -58 | 90 | 85 
MC5| -85 0.483 -60 | 90 | 93 


1.3.2 Stimulus Dynamics 


We show the results of our simulation with a stimulus for each model configuration 
(MC1-MC5) in Figures 1.3 and 1.4. These figures show the temporal and spatial 
changes, respectively, in ionic concentrations in the ICS and ECS. We also present 
the time evolution of the changes in the volume fractions of the ICS and ECS in 
Figure 1.5, subfigures a and b, respectively. 
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Fig. 1.3: Temporal dynamics in changes in ionic concentrations in ICS and ECS 
for different model configurations (MC1-MC5). The system is simulated for 60 s, 
with stimulation activated between [10 s, 20 s]. We measure the system at L = 
1.50-107^ m. The change in HCO; concentration is only shown for MC3 and MCS. 
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Fig. 1.4: Spatial profiles of the change in ionic concentrations in ICS and ECS for 
different model configurations (MC1-MC5). The system is simulated for 60 s, with 
stimulation activated between [10 s, 20 s]. We measure the system at the end of the 
stimulation, t = 20 s. The change in HCO; concentration is only shown for MC3 and 
MCS. 
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Fig. 1.5: Changes in volume fraction œ, as a function of time in ICS and ECS for all 
model configurations (MC1-MC5). 


We anticipate that the ionic concentrations in the compartments will remain 
constant until the stimulus is applied. This is seen in Figure 1.3 for all model 
configurations. For example, in subfigure a, we can observe that the intracellular Na* 
concentration remains the same until the stimulus begins in all model configurations. 
We can observe that the behaviour of MC3 and MCS are very similar throughout the 
temporal evolution of all ionic species (Figure 1.3, subfigures a - h). The increase 
of [K*]. and [K*]; can be seen most prominently in MC1 and MC2, and we see a 
similar effect of stimulation on these models in terms of extracellular and intracellular 
[Na+]. MC4 has the lowest magnitude change in [CI ]. among the models that do 
not simulate HCO; dynamics (MC1, MC2 and MC4). The model configurations 
that contain the NBC cotransporter for HCO; (MC3 and MC5) show a much lower 
influx of CI. 

In terms of spatial dynamics, we can see in Figure 1.4 that the behaviour of 
models MC3, MC4 and MCS is similar for all ions concentrations in ICS and ECS, 
except for extracellular Cl” concentration as shown in subfigure f. 

In both subfigures a and b, the volume fraction remains constant before applying 
the stimulus, as expected. After the stimulus was applied, we can see an increase in 
the volume fraction in the intracellular space and a decrease in the extracellular space. 
In fact, in response to ionic changes, the volume of the intracellular compartment 
increases by more than 6% for MC3 and MCS. The two model configurations where 
we can see the smallest increase are MC1 and MC2. Conversely, in response to 
the ionic changes, the volume of the extracellular compartment decreases to 12.5% 
for MC3 and MCS. As before, the model configurations where we see the smallest 
decrease are MC1 and MC2. 

The total intracellular and extracellular volume-weighted fluid velocities as ob- 
served at the end of stimulation (at t = 20 s) are shown in Figure 1.6, subfigures 
a and b, respectively. All fluid velocities are multiplied by the volume fraction of 
each compartment to obtain a volume-weighted result. The maximal ICS volume- 
weighted fluid velocities for all model configurations are also summarised in Table 
1.5. Here, we see that MC1 and MC2 have similarly low maximum water velocities 
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in ICS, while MC3, MC4 and MCS have similarly high maximum values (with MC3 
and MCS having the maximum ICS water velocity of 16 um/min). We also note 
that MC3 and MCS have almost identical fluid velocity plots. All models show fluid 
flowing in the same direction; positive intracellular fluid velocities for x > 150 um 


and negative intracellular fluid velocities for x « 150 um and vice versa in the case 
of extracellular velocities. 
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Fig. 1.6: Spatial profile of the total volume-weighted fluid velocity in the ICS and 
ECS for all model configurations (MC1-MC?5). 


Table 1.5: Maximal total compartmental volume-weighted fluid velocities for dif- 


ferent model configurations (MC1-MC5) during neuronal stimulus measured at the 
end of the stimulation (t = 20 s). 


Uj, max 


MCI 0.96 
MC2 1.1 
MC3 16 
MC4 14 
MC5 16 


We also present the osmotic and hydrostatic components of the ICS fluid velocity 
for MC1 and MC3 in Figure 1.7. When comparing the osmotic and hydrostatic 
components of the ICS fluid velocity between the model with the highest maximum 
velocity, MC3, and the model with the lowest maximum velocity, MC1, we see that 


the higher total ICS velocity in MC3 is due to the ratio of the osmotic contribution 
to the hydrostatic contribution. 
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Fig. 1.7: Comparing the volume-weighted osmotic and hydrostatic components of 
intracellular fluid velocity for MC1 and MC3. 


1.4 Discussion 


We have expanded on a previous computational model for predicting ionic electrod- 
iffusion and fluid dynamics in brain tissue by adding glial KCC1, NBC and NKCCI 
cotransporters. Further, we assess how the various cotransporters affect character- 
istics such as the energy required to maintain equilibrium, spatial and temporal 
dynamics of ions, and water velocities during neuronal activity. Our findings show 
that with no stimulus, MC5 has a lower Na*/K* pump flux value than the other 
models, from which we can infer that the astrocyte expends less energy with this 
configuration of membrane mechanisms. When the stimulus is applied, MC1, MC2 
and MC4 showed the most similarity for all the ion concentration dynamics. We also 
found that MC3 has the highest total intracellular of 16 tum/s. 


1.4.1 Calibration and Zero Stimulus Dynamics 


When there is no stimulus, we found very similar steady state ion concentration values 
across our model configurations, with the exception of intracellular and extracellular 
Cl” concentrations (with [Cl]. and [CI ]; for MC2 the biggest deviation: not 
significant, p > 0.05). The Na*/K* pump is the only ionic transport mechanism 
in our model that requires chemical energy to move ions against concentration 
gradients. Therefore, we were interested in seeing whether the flux through the 
Na*/K* pump (jpump) would have different values at equilibrium for the different 
configurations of the model. We see from Table 1.4 that MCS has a significantly 
(p < 0.05) lower pump flux value than the other models. As MCS has both the 
NBC and NKCC1 cotransporters, this low pump flux can be attributed to their net 
action in equilibrium. The flux values for MC1, MC3, and MCA are identical to three 
significant figures, while MCS has the lowest flux value. This is interesting as MC5 
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has both cotransporters, NBC and NKCC1, and needs the least energy to maintain 
equilibrium due to its low pump flux value. 


1.4.2 Stimulus Dynamics 


Focusing on the temporal Na*/K* dynamics in Figure 1.3, we can see that there is 
a decrease in the ionic concentration of Na* in the extracellular space due to the 
efflux of sodium ions due to the neuronal activities represented by the stimulus. The 
decrease in [Na*]; as shown in Figure a is due to the shrinkage of the extracellular 
space (and thereby the expansion of the astrocytic volume) resulting in a lower 
concentration of Na* (in spite of the influx of Na* ions (as [Na*]. > [Na*]; during 
stimulation). 

Neuronal stimulation results in an increase of [K*]; at a greater rate than the 
decrease in [Na*];. This creates an outward electrostatic gradient. To balance this 
gradient, we see an inward flux of CI ions in all the model configurations. This 
behaviour is different in MC3 and MCS, where we see an exponential decrease 
in [CI ]e (subfigure f) as in the other models followed by a linear increase for 
the remaining period of stimulation and an exponential increase at the end of the 
constant stimulus period above the baseline concentration before decaying back to its 
equilibrium. We believe that this is due to the dynamics of HCO; (subfigures h and 
g), which sees an exponential influx, in a sense playing the role of Cl ^ as seen by the 
other model. Therefore, since Cl” has an outward flux from the astrocytic membrane 
during stimulation, the intracellular concentration of HCO; continues to increase 
exponentially throughout stimulation. Analysing subfigures g and h, intracellular and 
extracellular HCO; concentrations are quite similar in both MC3 and MCS. This 
means that the NKCC1 cotransporter does not influence the spatial dynamics in the 
changes in HCO; concentrations. 

The magnitude of the osmotic and hydrostatic flows in the ICS for MC1 are nearly 
equal and in opposite directions, as seen in Figure 1.7. However, this is not the case 
for MC3, where the osmotic flow is more powerful than the hydrostatic, leading to 
a higher net water velocity. The osmotic water velocity depends on the immobile 
ions. As we see similarly high net fluid flows for MC3, MC4, and MC5 compared to 
MCI and MC2, we believe that the cotransporters involving Na* ions could be one 
of the common factors among the former group that causes this. The presence of 
HCO; ions and its effect on the valency and number of immobile ions could lead to 
increased osmotic flow. This could explain the higher total water velocities observed 
in MC3 and MCS, as demonstrated in Table 1.5. 

We opted to adjust the initial concentrations of the ions to obtain new steady states 
for our model configurations, while empirically setting the membrane conductance 
for each ion channel. This was a deliberate choice, and one could instead make 
the case for setting the concentrations empirically and adjusting the membrane 
conductivity to reach equilibrium (as done by Østby et al. 2009 [3]). 


1 Microscopic Fluid Flow in the Brain 17 


Although the comparison between each of our model configurations is useful, 
validation with empirical data could provide more biologically relevant information. 
While we are able to determine the effect of different mechanisms and their effect 
on shrinkage, a future goal is to compare these results with in vivo experiments and 
measurements of volume changes along with knockout experiments to identify the 
contribution of each membrane mechanism to determine the validity of our model. 
Future objectives to extend this model include incorporating additional membrane 
processes found in astrocytes, such as active water transport through the NKCCI 
channel and the inward rectifying K* current described in [2]. 


1.5 Conclusion 


In conclusion, our research underscores the role of astrocyte networks and glial 
cotransporters in the regulation of brain volume and waste clearance. Through our 
computational model within the Kirchhoff-Nernst-Planck electrodiffusive frame- 
work, we have identified the NBC cotransporter as a key player, exhibiting the 
highest fluid velocities. These findings highlight the implications of targeting glial 
cotransporters to modulate brain volume homeostasis and waste clearance, offering 
new directions for neuroscience research. 
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Chapter 2 S 


Computational Modeling of Ephaptic coupling 
in Myelinated and Unmyelinated Axon Bundles 
Using the EMI Framework 


Alessandro Gatti, Ramón Nartallo-Kaluarachchi, Abhinav Uppal, and Pietro 
Benedusi 


Abstract This report examines the Extracellular-Membrane-Intracellular (EMI) 
framework for modeling action potentials along 3D axons. We investigate the effect 
of myelination and the potential for ephaptic coupling in this model. Additionally, 
we assess the convergence and stability of a range of Runge-Kutta time-stepping al- 
gorithms on simple geometries with manufactured solutions. We first analyze single 
axons and the influence of myelin on the speed of action potentials. Then, we use a 
3D geometry of nine cylinders to represent an axonal bundle and study the induced 
potential in the central axon in both myelinated and unmyelinated cases. Finally, we 
discuss the biological implications of ephaptic coupling and the importance of 3D 
modeling for precise simulations of spiking neurons. 


2.1 Introduction 


Hodgkin and Huxley's 1952 work [1] was a pioneering effort to quantitatively 
describe the chemical processes that cause excitatory neurons to fire. Since then, a 
variety of alternative models, such as FitzHugh-Nagumo [2] and Hindmarsh-Rose 
[3] ODE models, have been developed to model excitable cells. The emergence 
of neural networks has further led to the modeling of neurons as single points in 
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space with only time-dependent dynamics. In this report, we explore two biological 
scenarios in which spatial and temporal dynamics are necessary to accurately model 
the propagation and induction of action potentials. Specifically, we investigate the 
effects of myelination [4] and ephaptic coupling in axonal bundles [5]. To do this, we 
use the recent EMI model to incorporate spatial dynamics into a partial differential 
equation (PDE) framework [6]. 


2.1.1 Myelination 


The myelin sheath is an extended and modified plasma membrane that is wrapped 
around the nerve axon in a spiral pattern [7]. It is derived from and is part of Schwann 
cells in the peripheral nervous system (PNS) and oligodendroglial cells in the central 
nervous system (CNS). Each myelin-generating cell provides myelin for only one 
segment of the axon. The nodes of Ranvier, which are short portions of the axon left 
uncovered, are essential for the functioning of myelin. This myelin sheath increases 
the resistance of the axonal membrane, lengthening its electrical space constant, and 
thus facilitating signal transmission along the axon. Additionally, myelin decreases 
the capacitance of the axonal membrane, so that less charge (in the form of Na*) 
is required to depolarize the cell. Both of these effects increase the speed of action 
potential propagation. 


2.1.2 Ephaptic Coupling 


The term ephapse was proposed in 1941 by Arvanitaki [8] to describe neural struc- 
tures coming into contact or in close proximity, without coupling via anatomically 
differentiated synapses. Where Arvanitaki had been studying an experimental prepa- 
ration with two giant axons (from Sepia officinalis) forming an ephapse by making 
the axons touch (for 5 mm in 4 to 5 cm long axons), the term ephaptic coupling 
has since come to describe short-range coupling between non-contiguous neuronal 
membranes [9], alternatively described as ‘electric field effects’ [10]. 

Similarly to how spiking individual neurons can give rise to extracellular action 
potentials, network-level activity can alter the local electric environment in nervous 
tissue. This spatiotemporal variation in extracellular potential and its gradient elec- 
tric field feed back into the same network, inducing ephaptic coupling [5]. This 
electric ephaptic coupling can be contrasted with ionic ephaptic coupling, where 
local changes in extracellular ion concentrations can alter Nernst potentials [11]. 

Although ephaptic coupling between neurons can be considered a weak effect 
compared to chemical or electrical synapses (measurements of endogenous electric 
fields are on the order of a few mV / mm, see [12] for a summary of experimental 
data), it can still have network-level implications in healthy and pathological nervous 
tissue [10]. For instance, weak electric fields can entrain slow neocortical oscillations 
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[13]. Additionally, mathematical models of ephaptic coupling in axon bundles show 
induced synchronization of firing activity [14] and modulation of transmission delays 
[15]. 

In pathological scenarios, ephaptic entrainment has been implicated in neurode- 
generative disease models considering damaged neurons [16]. Ephaptic coupling 
is also hypothesized to play a role in pathologies resulting from demyelination of 
cranial nerves under compression, such as trigeminal neuralgia (facial pain) through 
the fifth cranial nerve [17] and vestibular paroxysmia (short episodic vertigo) from 
the eighth cranial nerve ([18], [19]). 

Thus, studying ephaptic coupling can improve our understanding of network-level 
feedback under healthy and pathological conditions. Modeling studies of ephaptic 
coupling also offer the translational benefit of providing an evaluation framework in 
the design of exogenous (transcranial) stimulation protocols [12]. Given the key role 
of extracellular potentials in ephaptic coupling, here we extend previous work on 
modeling axon bundles using the EMI framework [11] by contrasting unmyelinated 
and myelinated axons in a bundle. 


2.2 Methods 


2.2.1 The EMI Model 


The EMI model is a PDE framework that has recently been developed to simulate 
excitable cells, such as neurons and cardiac cells, from the first principles [6]. 
It divides the extracellular space, the cell membrane, and the intracellular space 
into distinct components. When both spatial and temporal dynamics are taken into 
account, it is suitable for modeling the effects of myelination and ephaptic coupling. 


For a single cell, denoted Q;, surrounded by an extracellular domain, Qe, the EMI 
model is given by the following coupled PDE-ODE, 


Ve ojVuj = f, in Q;, 
V -oceVue = g, in Qe, 
OeVue : ne = —0; Vu; - ni = In, at T, (2.1) 
V =U; — ue, at T, 
Ov 1 
a oS oF In = lion , t T, 
Bt Ca! ) i 


where ui, ue, and v are intracellular, extracellular, and membrane potentials, respec- 
tively, which are commonly given in mV. Furthermore, c; and c; are intracellular 
and extracellular conductances, respectively (typically in mS/cm), C,, is the mem- 
brane capacitance (typically in uF/cm?), and T denotes the cell membrane. n; and ne 
represent the outward-pointing normal vectors. The ionic currents through channels, 
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pumps, and exchangers at the membrane are denoted by Tjon and are typically given 
in 4A/cm?. A schematic representation of the model domain is given in Fig. 2.1. 


We assumed that the external boundary 0Q, is insulated, which leads to the following 
Neumann boundary condition 


Oe Vue : ne = 0, at Qe. (2.2) 


Ne QQ. 


Fig. 2.1: Illustration of an EMI model domain consisting of an extracellular domain, 
Qe, a cell membrane, I, and an intracellular domain, Q;. Ne, n; represent normal 
vectors and 0Q, represents the boundary of the extracellular space. 


2.2.2. Numerical Methods 


2.2.2.1 Multi-Dimensional Primal Formulation of the EMI Model 
In order to numerically integrate the EMI model, we derive a weak formulation [6]. 


In particular, we use the multi-dimensional primal formulation which is expressed 
as follows: find u; € V; = H! (Qi), ue € Ve = H! (Qe), Im € Q* = HV? (T) such 


that 
Í. ou Wid | vids. | fvidx, 
Qi r Qi 
7. ou, Wed - | Inveds = | gvedx, (2.3) 
Qe r Qi; 


[eins f winds- | nne nts [ hinds, 
r r pF r 


for all the test functions v; € Vj, Ve € V, and jm € Q* [6]. In this formulation, the 
time-dependent equations are discretized according to an implicit Euler scheme. The 
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known right-hand side h combines the previous transmembrane potential solution, 
vo, and the evaluation of the ionic current, J;5,, into h = vg — Co Atl isi: However, 
as discussed in the following section, this final equation can be modified according 
to any time discretization. 


2.2.2.2 Runge-Kutta Time-Stepping in the EMI Model 


As discussed above, the derivation of the weak formulation requires a choice of 
time discretization. Depending on the choice of the time-stepping algorithm, the 
final equation will differ. We consider the well-studied Runge-Kutta (RK) methods, 
a family of classical ODE integrators [20]. For a homogeneous differential equation 
of the form 


dy 


gro £00 Q4) 


an s-order RK scheme is a one-step method of the following form, 


Yni = Yn +h X biki, (2.5) 
i=l 
where, 
ki = fn). (2.6) 
ks = fs * h >) auk). (2.7) 
j=l 


Thus, the method is uniquely defined by the matrix vector pair A, b, where A = (aij) 
and b = (bj, ..., bs). We can also define the vector k = (k1, ..., k2). To make this 
method consistent with the order s, certain conditions on A, b must be satisfied, 
limiting the number of RK methods for each given order. When A is lower triangular 
with zero diagonal, the method is explicit; otherwise, it is implicit. We investigate 
the form of final EMI equation for a generalized RK scheme for both the passive and 
Hodgkin-Huxley conductance dynamics. 


2.2.2.3 The Passive Conductance EMI Model 


We consider the time-dependent equation from the EMI model with passive conduc- 
tance, corresponding to Ijo,(v) = v [6], 
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dv 1 
T = Cs (Im — v(t)). (2.8) 


Considering an RK method of order s, we discretise this equation as, 


Viel = Vet 33 biki, (2.9) 
i=l 
=v, +hb'k. (2.10) 
We further notice that 
pee aa (2.11) 
Cy = 


which we write in vector notation as, 


1 
k = cos - vi)u - hAk, (2.12) 


where u is a s-vector of 1s. Defining B = S and I to be the s-dimensional identity 
matrix, we assume I + BA to be invertible and write k explicitly as 


k= ec Un = v) BAY, (2.13) 


As such, our RK scheme can be written as 
vii = (1 - Bb (1+ BA) F'u)v, + Bb" (1+ BA) lly. (2.14) 


According to the weak formulation of the EMI model, this yields conservation law, 


fu ds- fue ds- B (I +BA)'u [In ds 
T T p 


- (1 - B8 Qe pAy'u) f v ds. (2.15) 
T 


Given the Butcher tableau, (A, b), one can use this formulation to calculate the 
function, 


eux (B) = Bb (I+ BA) !u, (2.16) 


which can then be used to implement any RK scheme on the EMI model. Further- 
more, we notice that grg (£) is related to the well-studied stability function of the 
RK family, which is defined [20], 


R(z) 21*zb'(I zAy !u, (2.17) 
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namely, the relation is, 


8RK(B) = 1 — R(-), (2.18) 


yielding conservation law, 


fu ds- fue ds (RC) e ) [In as= RCD [ v ds. (2.19) 


The stability function can be written as 


det(I — zA + zub") 
det(I — zA) 


R(z) = ; (2.20) 
where for explicit scheme, due to their lower triangularity, the denominator is 1. 
Therefore, for explicit schemes, the stability function R and therefore ggg is a 
polynomial of degree s, which we denote p; (5) that has form, 


ps(B) = X CD" g"G(A, b, n), (221) 
n-l 
with, 
s ig-l ik-ı-1 in-2-1 
GU bn) S bg >) ase D, ee > dues (22D) 
ig=n ij-n-l iyzn-k in-j=n-1 


for n = 2, ..., s and extended to n = 1 by, 
@(A,b, 1) = > bi, (2.23) 
ig=1 


This yields conservation law, 


ji i die J ue ds — ps(B) J Im ds = (1— ps(B)) i vids. (2.24) 


Following this process, we can derive the functions ggk associated with a range of 
explicit and implicit RK schemes. In addition, we note that for explicit RK schemes 
of order p, we have the following approximation, 


g^. oo zP 
R(z)=1l+zt—ot—4..+¢ — + O(z?*4), (2.25) 
2! 3! p! 


therefore R(z) = e7 + O(zP*!). Using this property, we can define the following 
scheme that approximates any explicit, p-order RK scheme to order p, 
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Table 2.1: Common Runge-Kutta schemes and the associated functions grx (£). 


Scheme Abbrev. Order Implicit/ gnk(B) 
Explicit 
Explicit Euler EE 1 E B 
Explicit Midpoint EMP 2 E g- g 
Runge-Kuta4 ^ RK4 4 E dicÉalgted 
Implicit Euler IE 1 I 5 
Trapezoidal Rule TPR 2 I $E 
Exponential function EF peN E 1- e£ 


fu ds- fue ds + (e P + ) fm as=e f v ds, (2.26) 
F r T r 


which is referred to in Table 2.1 as “Exponential function” (EF). This scheme 
represents the exact solution of the linear time-dependent equation. 

As shown above, using both higher-order and implicit schemes in the passive 
conductance model can be done via direct calculation, and therefore at no additional 
computational cost. 


2.2.2.4 Numerical Convergence of RK Schemes 


In order to compare the convergence and stability of the different RK schemes, we 
consider the simple square unit domain shown in Fig. 2.2. 


(0, 1) (1, 1) 


(0, 0) (1,0) 


Fig. 2.2: Domain of integration. Qe, O;, [' represent the extracellular space, the 
intracellular space and the membrane, respectively. 


Using the manufactured solution method [21], we define the source functions 
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f= f@ y,t) = -87° sin(2zx) sin(2ay)(1 +e), (2.27) 
g = g(x, y, t) = -87° sin(2zx) sin(2zy). (2.28) 


In this case, the exact solution is given by, 


uj(x, y, t) = (1 e *) sin(2zx) sin(2zy), (2.29) 
uc (x, y, t) = sin(2zx) sin(2zy), (2.30) 
v(x, y, t) = e ' sin(2zx) sin(2zy). (2.31) 


We calculate the L? error between the numerical solutions, aN AE a E pyar using 


the finite element method (FEM) and the exact solutions above. We present the two 
error terms, 


"N.N NU 
e^t — Alu; — aA + Jue — ^p, (2.32) 


gh At = ilv _ pNAt||2, (2.33) 


We assess the convergence with respect to both temporal and spatial resolution. 
For the spatial convergence test, we fix At = 0.01/64 and vary the mesh resolution 
parameter, N = 16, 32, ..., 256, which equals the number of intervals in each direction 
of spatial discretisation of the domain, as shown in Fig. 2.5. For the temporal 
resolution, we fix N = 512 and perform two experiments. Firstly, we use At = 0.01/n 
for n = 0.1,0.5,1, 2, 4, 8, 16, 32, as shown in Fig. 2.6. 


2.2.2.5 The Active Conductance Model 


For the EMI model with active conductance, the time-dependent equation becomes 
[6], 

dv 1 

dt C. (Im — Tion(v(t))) , (2.34) 
where Jion(v(f)) is a non-linear function, for example of Hodgkin-Huxley type. In 
this case, explicit expressions cannot be derived when applying RK schemes. Given 
an implicit RK scheme with Butcher tableau, (A, b), the scheme is again, 


Vai Span bi biki, (2.35) 
i=l 
=v, thb'k. (2.36) 


However, the equations for k; yield a non-linear system of equations to be solved at 
each time-step, 
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1 S 
ki = c, Un = Tion(Vt — ee) (2.37) 
which we write in vector notation as, 
1 
k = — (Im — Fion(vr, A, k)), (2.38) 


Cm 
where, 
Tion(vt + ej Ak) 
Fion(vi, A, k) = : ^ (2.39) 
Tion(v; + e1 AK) 


where (e;); , are the standard basis vectors in R*. 


Most frequently, one considers an active conductance of Hodgkin-Huxley type [1], 


Tion(V) = gua(v) (v — Ena) + gk (v) (v — Ex) + gL (v) v — EL) + g(v)(v - Ena), 
(2.40) 


where the gating variables, gion, for each ion, are non-linear functions of v given by 
the solution to an ODE. This form allows us to define semi-implicit schemes for this 
model. At each time step, the gating variables are updated according to a numerical 
method, typically the Rush-Larsen method [22]. To define a semi-implicit scheme, 
we update the gating variables with membrane potential v; rather than the value 
vt+h 2-34 a;jkj. This now yields a simplified system of equations for k;, 


S 


1 S 
k; = c, Un — gua(vi) (v; + ELI — Ena) - gk(vi) (vi + haut — Ek) 
- a.) +h X agk; - E) - 8(v:)(vr +h J agk;- Ex) — (2.41) 
j=l j=l 


which is linear. In vector notation, the system becomes 


1 
k= T [m — lios(vi)]u — BGAK, (2.42) 


where G = gua(v;) + gk(vi) + gi (vi) + g(v;). This has solution, 


k= e Un - li (ve) | + BGA) !u. put 


Thus, the scheme becomes as follows: 
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Viet = Vr Bb! [Im — lis (vi))] (I + BGA) 'u. (2.44) 


Under the weak formulation this gives us conservation law, 


fu ds- fue ds- Bb (1 +BGA)'u [ In ds 
T r r 


= fv- ion (vi) b (1 + BGA) u ds. (2.45) 
Dr 


However, this inverse matrix that must be calculated equates to solving a linear 
system at each time step. 


2.2.3 Model Setup and Geometries 


We employ the multi-dimensional EMI formulation (eq. 2.1) with an implicit Euler 
scheme (eq. 2.4). We use a Hodgkin-Huxley active conductance model [1] integrated 
with a Rush-Larsen scheme [22]. The source terms f and g in eq. 2.1 are set to 0 [6]. 
The model parameters were chosen as shown in Table 2.2. The conductance values 
for the unmyelinated axons correspond to the giant squid axon [1]. For the myelinated 
axons, we scaled the corresponding unmyelinated conductances such that the total 
number of ion channels would match for the two axon types. The axon length, radii, 
and node lengths were chosen to make the geometry and results easier to visualize. 
The synaptic conductance g(v;) in Iion represents a stimulation current applied to a 
synapse of length node on one edge of unmyelinated and myelinated axons. 


2.2.3.1 Single Axon with and without Myelination 
Fig. 2.3 illustrates the geometries we employed in our simulation to compare the ve- 


locity ofthe action potential in a single axon for the unmyelinated (a) and myelinated 
(b) scenarios. Table 2.2 contains the parameters used for the axons. 
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Fig. 2.3: ParaView renders of the gmsh geometries for (top left) unmyelinated, and 
(top right) myelinated single axons, where (bottom) shows a schematic representation 
of both axon types using parameters from Table 2.2. 


Table 2.2: Parameters used for single axon simulations with and without myelin. 


Parameter 


Axon length L (mm) 
Intracellular radius rj, (mm) 
Extracellular radius rex (mm) 
Myelin thickness tye) (mm) 
Number of nodes of Ranvier 
Length of node Inode (mm) 
Internodal distance /,,, 4; (mm) 
gr (mS/cm?) 

&Na(mS/c m^) 

gk (mS/cm?) 

Er (mV) 

Ewa (mV) 

Ek (mV) 


Value 
Unmyelinated Myelinated 
10 10 
0.2 0.2 
1 1 
- 0.2 
- 10 
- 0.1 
- 1 
0.3 3 
120 1200 
36 360 
-54.38 -54.38 
54.8 54.8 
-88.98 -88.98 


2.2.3.2. Axon Bundle 


For our simulations considering bundles of unmyelinated axons, we arranged nine 
axons in a 3 x 3 grid. A cross-section showing this spatial arrangement is depicted 
in Fig. 2.4. Each axon in the bundle was matched to use the same parameters and 
geometries as the single axon simulations described previously. The same arrange- 
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ment and parameter matching to single axon simulations was also applied to a 3 x 3 
bundle of myelinated axons. 

For both the unmyelinated and myelinated axon bundles, the eight axons on the 
periphery were stimulated using synapses of length /node (as in the single axon 
simulations), while the central axon in the bundle was not stimulated externally. 
Membrane voltage changes in this central axon were examined to look for evidence 
of ephatic coupling. 


Fig. 2.4: Arrangement of nine axons forming a bundle. Q; denotes the intracellular 
space of each axon, Qe denotes the extracellular space, and I denotes the membranes. 


2.3 Results 


2.3.1 Convergence 


10? 
> 
[] 
5 10? 
o 
N 
Ll 
1074 
10? 10? 
N: Spatial mesh resolution parameter N: Spatial mesh resolution parameter 


Fig. 2.5: Convergence with respect to mesh size with fixed time-step. Left: e,, Right: 
ey. 
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Fig. 2.6: Convergence with respect to temporal resolution for small time-steps and 
fixed mesh size. Left: e,, Right: ey. 


Firstly, we consider the convergence results for the square-in-square, passive conduc- 
tance, manufactured solution discussed previously. We integrate the system using 
FEniCS and the RK schemes listed in Table 2.1. In Fig. 2.5, we display the error for 
a sequence of increasing spatial (mesh) resolutions for both the intra-extracellular 
potential and the membrane potential. As shown, all the methods converge com- 
parably with respect to mesh resolution, as they only modulate the time-stepping 
error. Next, we consider a series of decreasing time steps. As can be seen in Fig. 2.6, 
higher-order methods are limited by spatial discretisation, but all methods appear 
to converge to the exact solution. It is clear that the higher-order methods, EMP, 
RK4 and TPR have lower errors than the first-order methods EE and IE. For these 
higher-order methods, in order to assess their convergence and their stability, we 
find the time steps at which they obtain similar error values. By increasing the time 
steps so dramatically, we can limit the effect of the spatial discretisation error and 
compare the methods more effectively. As shown in Table 2.3, the least stable of 


Table 2.3: Errors for large time-steps for higher order RK schemes. 


At RK4 At EMP At TPR 
eu ey eu ey eu ey 


2.5 1.03x 1075 5.62x107 1.6 1.03x107? 2.15x 1074 3.0 1.60x 10? 2.50x 1074 
2.6 2.99x107> 4.70x107^ 1.8 1.64x107? 2.15x 1074 3.2 1.30x 107^ 2.17 x 10? 
2.7  1.74x107? 2.90x 107? 2.0 1.65x 107! 2.78 34 112x102 1.87x107 


the higher order methods is EMP, followed by RK4 and finally TPR, the implicit 
scheme. Importantly, the higher order explicit schemes had much lower error yet are 
less stable. This means that, depending on the constraints of the particular problem, 
one should opt for an implicit scheme for greater stability or an explicit scheme for 
greater accuracy. 
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2.3.2 Speed of Action Potential Propagation in Single Axons 


We conducted a simulation to compare the propagation speed of an action potential 
in a myelinated axon with that of an unmyelinated one. Fig. 2.7 shows the membrane 
potential in both cases recorded at a single point on the membrane and at the end of 
the axon. The action potential profile is similar in both cases. The key difference is 
that the potential propagates faster in the myelinated axon; therefore, it arrives at the 
end of the axon earlier, as shown in Fig. 2.7. 


] 
—— Unmyelinated 


40 —— Myelinated 


Potential (mV) 


Time (ms) 


Fig. 2.7: Potential of myelinated and unmyelinated axons over time measured at the 
node located at the end of the axon membrane. 


Our simulation supports the biological theory that argues that myelination causes 
faster propagation of action potential [23]. This further validates this approach to 
modeling myelination by modulating the conductance values in the nodes of Ranvier. 


2.3.3 Ephaptic Coupling in Unmyelinated and Myelinated Axon 
Bundles 


Next, we simulate the 9-cylinder axonal bundle with myelinated and unmyelinated 
axons, while stimulating only the 8 peripheral axons. We then measured the induced 
activity in the unstimulated central axon. 
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Potential (mV) 


o 2 4 6 8 10 0 2 4 6 8 10 
Time (ms) Time (ms) 


(a) Comparison between the action potential (b) Comparison between the action poten- 
of the central and peripheral axon in the un- tial of the central and peripheral axon in the 
myelinated case. myelinated case. 


— Myelinated 
— unmyelinated 


Potential (mV) 


Time (ms) 


(c) Comparison between the action poten- 
tial of central axons in the unmyelinated and 
myelinated cases. 


Fig. 2.8: Plots showing the action potentials induced in the central axon alongside (a 
chosen) peripheral axon in both the myelinated and unmyelinated cases. 


In Fig. 2.8 (subfigure 2.8a), we show that in the unmyelinated case, ephaptic coupling 
did occur as an action potential was induced in the central axon. There is a noticeable 
time-lag between the spike in the peripheral and central axons. In the myelinated 
case, shown in Fig. 2.8 (subfigure 2.8b), ephaptic coupling also occurred with an 
induced potential in the central axon. There is no such time-lag in the case of the 
myelinated axons. The different profiles of the induced action potentials in each 
case are highlighted in Fig. 2.8 (subfigure 2.8c), which shows the unmyelinated and 
myelinated central axons side by side. 
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Fig. 2.9: Action potentials induced in the central axon of the unmyelinated bundle 
for different time-stamps. Arrows depict the splitting of the induced potential into 
forward and backward propagating impulses. 


Interestingly, in our simulations we have observed a backpropagation phenomenon 
in the central unmyelinated axon as one can see from Fig. 2.9. 


Our simulations show that ephaptic coupling is possible in the EMI model, for both 
myelinated and unmyelinated axons. Furthermore, we have shown that myelination 
severely increases the speed of the propagation of action potentials under direct 
stimulation and indirect induction. 


2.4 Discussion 


In this report, we presented a computational approach to modeling myelination and 
ephaptic coupling of axons using the EMI framework [6]. As part of this analy- 
sis, we studied a range of time-stepping algorithms for both the passive and the 
Hodgkin-Huxley conductance models, obtaining an analytical expression for the 
weak formulation for generalized RK schemes. As a result, we can integrate the EMI 
model with any RK scheme. We concluded that higher-order schemes are advan- 
tageous in terms of both stability and accuracy, with explicit methods being more 
accurate and implicit methods being more stable. Our computational simulations 
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also validated the hypothesis that myelination speeds up the propagation of action 
potentials. Finally, we showed that ephaptic coupling was possible in the EMI model 
for both myelinated and unmyelinated axons, where myelination accelerated the 
propagation of the induced action potential. 

Our work further supports the notion that biophysical models of excitable cells 
should take into account spatial geometry and dynamics to model more specialized 
biological phenomena. Both myelination and ephaptic coupling are important fea- 
tures of neuronal function and show pathology in disease, such as myelin loss in 
multiple sclerosis [24] or reduced ephaptic entrainment in neurodegenerative condi- 
tions [16]. For this reason, our report represents an important step towards complete 
models of excitable cells and the phenomena that can arise. 

Although the EMI framework is a detailed and biophysically realistic model, 
further work using the more detailed KNP-EMI model [11], would allow for a more 
complete analysis of the parameters that modulate the emergence of induced action 
potentials, such as the distance between the axons in the bundle. Parameter values 
as measured in demyelinated axons could also be considered for modeling ephaptic 
coupling under pathological conditions [25]. 

Furthermore, this report is a proof-of-concept showing that modeling ephaptic 
coupling and myelination is possible within this model, yet does not claim to replicate 
the exact dynamics seen in biological experiments. Further studies could confirm 
that the model outputs are consistent with theoretical predictions, that is, from the 
cable equation [26], or from experimental results. 
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Augmentation of Cardiac Ischemic Geometry for 
Improving Machine Learning Performance in 
Arrhythmic Risk Stratification 
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Abstract Ventricular arrhythmias frequently occur as a complication of myocardial 
infarction (MI), due to significant changes in the heart's structure and electrophysi- 
ology. If left untreated, these alterations may lead to sudden cardiac death (SCD). It 
is therefore critical to evaluate risk prediction accurately in post-infarction patients 
to enable early intervention and improve patient outcomes. This work introduces a 
novel approach to improve arrhythmia risk assessment in post-infarction patients. 
We propose a new pipeline to build physiologically realistic image-based models 
of patient hearts, producing more realistic meshes compared to publicly available 
pipelines. We generate a library of 90 cardiac geometries of MI patients and use 
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these cardiac models to estimate likelihood of reentry using electrophysiological 
(EP) simulations. However, due to the computationally expensive nature of this ap- 
proach, we also introduce a data augmentation pipeline to train a machine learning 
(ML) model for risk stratification, enabling accurate and real-time prediction of the 
simulation outcomes. Our trained ML model achieved an accuracy of 88.0% and F1 
score of 48%, with a prediction time of 0.01 seconds per case (compare with approx- 
imately 5 hours per case for EP simulations). In conclusion, the work presented here 
improved the accuracy of personalised biventricular geometries, introduced a novel 
data augmentation approach for scar distribution, and decreased prediction time of 
risk of arrhythmias post-MI by more than five orders of magnitude. 


3.1 Introduction 


Ventricular arrhythmias are a common consequence of MI and are associated with a 
significant increase in SCD in post-MI patients [1, 2]. Infarction occurs when blood 
flowing to a portion of the myocardium decreases or stops flowing altogether, typi- 
cally due to a narrowed or completely occluded coronary artery. A rapid reparative 
process is initiated to rebuild the infarcted myocardium and preserve the structural 
integrity of the ventricles, resulting in the formation of a scar [3]. Ventricular scar 
tissue is characterised by dense fibrosis, comprising of collagen and fibrocytes, as 
well as regions containing surviving myocyte bundles. The heterogeneous presence 
of dense fibrosis can lead to conduction block, delineating the borders of potential 
reentry circuits. Fibrosis also contributes to the development of slow tissue conduc- 
tion, a pivotal factor in the formation of reentry mechanisms [4]. The combination of 
these factors creates a vulnerable environment, acting as a substrate for ventricular 
arrhythmias. The role of fibrotic tissue in initiating and maintaining ventricular ar- 
rhythmias is indeed complex and remains an active area of research [5, 6, 7]. Timely 
and accurate prediction of these arrhythmias is crucial to enable early intervention 
and improve post-MI patient outcomes, and as such represents an important clinical 
challenge to address. 

The occurrence of arrhythmic events in patients at risk may be prevented by 
an implantable cardioverter-defibrillator (ICD), making these devices the primary 
treatment of choice for high-risk post-MI patients [8]. However, ICD implantation is 
costly and associated with adverse effects [9]. Therefore, there is value in developing 
non-invasive methods for prospective identification of post-MI patient subgroups 
with the highest risk of arrhythmia. 

Image-derived mechanistic computational models of cardiac electrical function 
have proven successful in predicting arrhythmic risk [10, 11]. Despite their effec- 
tiveness, these methods rely on complex, computationally expensive models. As a 
result, there is a growing need for more efficient approaches. By replacing mech- 
anistic simulations with a statistical model, it is possible to learn the relationship 
between the inputs and outputs of a physical process, based solely on data and with 
no phenomenological information. Deep learning (DL) models have emerged as a 
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promising solution, due to their ability to analyse large amounts of data and identify 
complex patterns [12, 13, 14, 15]. While these statistical methods reduce computa- 
tional cost compared to mechanistic simulations, achieving accurate results requires 
a substantial amount of training data [16]. 

Data augmentation can help address this issue by increasing the size and com- 
plexity of the training dataset, thanks to manipulations performed on the existing 
data. Augmenting the dataset enhances the model's ability to generalise to unseen 
cases, by exposing it to more heterogeneous cases and mitigating overfitting [17]. 
To this end, biophysical simulations of cardiac electrical function can be used to 
generate large amounts of data within physiological expectations. Thus, DL models 
are able to better capture process intricacies at a local scale and better represent 
population heterogeneities at a global scale, improving their accuracy in performing 
personalized predictions of arrhythmic risk. 

This study aims to enhance risk stratification of post-MI patients using augmented, 
biophysically-accurate image-based models of the heart. Building on the work by 
[14], we propose a novel meshing pipeline to obtain realistic anatomical models based 
on patient magnetic-resonance (MR) images, including ventricular scar location and 
fiber orientation. To increase the variability of our population, we augment the 
original anatomical models obtained from clinical MR images by changing scar 
characteristics including size, thickness and location. We then perform personalised 
electrophysiological simulations on over 100 different modelled cases to obtain 
activation sequences in the ventricles using varying conduction velocities to reflect 
tissue properties in healthy tissue and scar. Subsequently, we determine the risk of 
arrhythmia of each case based on the absence or presence of reentry observed in the 
simulation results. 

Understanding the pathophysiological link between anatomical properties of the 
healthy myocardium and scar and ventricular arrhythmias can inform approaches 
to identify patients at risk. Using the augmented dataset and simulation outputs, 
we train a graph convolutional neural network (GCN) to predict the occurrence of 
reentry based on regional properties of the myocardium. Our model is able to emulate 
the simulation results within 88.0% accuracy and improve prediction time by five 
orders of magnitude, thus enabling rapid, personalised risk stratification. 


3.2 Methods 


3.2.1 Dataset 


The available data consists of standard short axis late gadolinium enhanced magnetic 
resonance images (LGE MRI) provided by Rigshospitalet in Copenhagen, DK, which 
counts 30 confirmed cases suffering from first-time MI. The reader can refer to [14] 
for further details. 
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3.2.2 Data Augmentation 


3.2.2.1 Ventricular Mesh Generation 


For each scan, the LGE MRI stack was imported into ADAS 3D (version 2.12.0) 
and the slices were automatically aligned using their pre-trained model. The left 
ventricular (LV) epicardium and endocardium contours were automatically generated 
using the ADAS 3D auto-segmentation model with only minor manual adjustments 
needed. For the right ventricle (RV), the epicardium and endocardium were manually 
segmented and the contours were then generated. 


— LV epicardii 

_ LV endocardium 
— RV epicardium 
+++ RV endocardium 


Fig. 3.1: Left: Short axis LGE MRI. Right: Contours of the epicardium and endo- 
cardium of the left ventricle (LV) and the right ventricle (RV). 


The extracted contours (Figure 3.2A) were refined using an automated pre- 
processing script based on Blender (version 3.01). The refinement ensures that the 
ventricular base is consistently flat (Figure 3.2B). The resulting surfaces (ventricular 
base, LV epi- and endocardium, RV epi- and endocardium and septum) represent a 
single closed surface of a biventricular mesh (Figure 3.2C). 
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Fig. 3.2: (A) Surface mesh generated from epicardial and endocardial contours of 
the left and right ventricles (green) and surface mesh of the scar generated from 50% 
FWHM (red). (B) Ventricular mesh with flat base. (C) Closed-surface biventricular 
mesh and smoothed scar mesh. 


The surfaces are combined, scaled and rendered as a single 3D tetraheadral mesh 
using gmsh (version 4.8.4) [18]. Based on considerations of [14] we chose a spatial 
resolution of 400 pm for each generated mesh. At this point, the generated 3D mesh 
does not include any cardiac muscle fiber information. To generate fiber information 
for a given mesh an implementation based on the works of [19] was used. Finally, 
in order to be able to perform openCARP simulations [20], the resulting mesh is 
converted to the respective CARP format using meshtool [21]. 


3.2.2.2 Scar Mesh Generation 


After generating the LV epicardial and endocardial contours in ADAS 3D, the LV 
scar tissue was characterized using full width at half maximum (FWHM) threshold 
method. For every LGE MRI scan, at least 3 distinct scar geometries were generated 
using different FWHM percentages: 30%, 50% (Figure 3.2A), and an arbitrary 
percentage value between 20% and 70%, other than 30% or 50%. Subsequently, 
each scar surface mesh underwent refinement in Blender to ensure that the scar 
tissue is contained within the left ventricle and that small scar chunks with a surface 
smallerthan 10 mm? are discarded. Furthermore, the refinement process improves the 
quality ofthe scar surface by eliminating any artifacts and smoothing the scar surfaces 
(Figure 3.2C). Finally, the scar surfaces are incorporated into the 3D tetrahedral mesh 
generation process. 

In order to extend the scar population per mesh, an artificial scar generation algo- 
rithm was developed which injects scar tissue into given biventricular geometries. 
The algorithm solves an Eikonal wave equation for a given set of either manually 
or randomly selected source points within the mesh [22]. The volumes contained 
within the propagated wavefronts, after an adjustable period of time, are labeled as 
new scar tissue (Figure 3.3). Furthermore, inner regions of generated scar tissue are 
labeled as infarct tissue, while the outer regions are labeled as border zone. The 
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velocity of the wave propagation can be controlled in order to generate complex scar 
shapes. 


Fig. 3.3: Artificial scar generation illustration with example. Based on source points 
(yellow) a propagating wave is computed. The enclosed volume represents new scar 
tissue with border zone (light red) and ischemic tissue (dark red). 


3.2.2.3 Mesh Feature Extraction 


After generating the LV epicardial and endocardial contours in ADAS 3D, the 17 
American Heart Association (AHA) segments were computed for the LV. The AHA 
segments are a widely accepted clinical convention that define different anatomical 
regions of the LV [23]. For each segment, we calculated the segment mass, LV 
wall thickness, scar mass within the segment, and scar percentage. In addition to 
the segments, the region between the LV endocardium and epicardium was divided 
into 10 layers (with Layer 1 corresponding to the endocardial surface and Layer 10 
corresponding to the epicardial surface). For each segment, we also calculated the 
layer surface area, scar area, and scar percentage for each of the 10 layers within the 
segment. 


Layer 10 
mE. Se 


3 


Healthy 


Fig. 3.4: Left ventricle anatomy subdivided into 17 AHA segments and 10 layers 
from endocardium to epicardium. Scar threshold at 50% FWHM. 
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3.2.3 Electrophysiological Simulations 


Building upon the work of [14], we incorporated the Ten Tusscher model [24] to 
simulate electrical wave propagation in healthy cardiac tissue. By manipulating ionic 
conductances and tissue conductivity reduction along transverse and longitudinal 
axes, our model captured the intricate progression of damage occurring post-MI. We 
employed the same parameters for healthy tissue and ischemic regions as those used 
in [14]’s work. 


3.2.3.1 Stimulus Point Selection 


For the selection of the stimulation points, we follow a different strategy than [14]. 
Instead of selecting a stimulation point per AHA segment, we decided to select 
stimulation points which are positioned near the scar. [25] showed that premature 
beats in acute ischemia first occur in healthy myocardium adjacent to ischemic 
regions; stimulation points near the scar tissue are therefore more likely to induce 
reentry. Furthermore, the selection of stimulation points follows a similar strategy as 
the artificial scar generation approach presented in Section 3.2.2.2. Specifically, the 
border zone of each scar is used as a starting point for the Eikonal wave propagation 
algorithm. The volume enclosed within the propagated wavefront after a short period 
of time represents tissue which is close to the scar (Figure 3.5). Within this volume, 
we select a number of points and ensure that they are spaced apart by a minimum 
distance. 


Fig. 3.5: Stimulation point selection illustration with example. Starting from the 
border zone regions (light red) a propagating wave is computed. The enclosed 
volume (light yellow) represents a region near the scar tissue. Within this region up 
to 17 stimulation points are selected. 


3.2.3.2 Vulnerability Simulation Protocol 


We conducted the simulation protocol at the different locations (see Section 3.2.3.1) 
pointed out in [14]. In summary, each simulation involved the delivery of pacing 
stimuli, known as S1. Initially, these S1 stimuli were administered with a fixed cycle 
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length of 350 milliseconds, followed by an S2 stimulus 200 milliseconds later. If 
no irregular heart rhythms were observed during this phase, we gradually reduced 
the time interval between S1 and S2 in 10-millisecond increments. We continued 
this process until we either detected abnormal reentrant circuits within the heart, or 
until the S2 stimulus failed to propagate. In cases where the S2 stimulus failed to 
propagate, we introduced an S3 stimulus 250 milliseconds after the last successful 
S2 delivery, and followed the same procedure as before. If reentry was not detected 
even after the S3 stimulus, we introduced an additional S4 stimulus 250 milliseconds 
later, following the same protocol as with the S2 and S3 stimuli. Finally, the resultant 
outcomes were categorized as follows: absence of reentry, temporary reentry, or 
persistent reentry. 


3.2.4 Arrhythmic Risk Prediction Model 


The electrophysiological simulations described above provide an accurate estimation 
of ventricular activation times and, if present, re-entrant circuits in each subject. Sim- 
ulation outputs are obtained using a fully phenomenological, multi-scale modelling 
pipeline that requires high performance computing (HPC) resources. As such, this 
approach is computationally expensive. From a clinical perspective, this limitation 
creates a need for faster methods that can reliably estimate of the outcome of interest 
whilst avoiding the requirement of HPC resources, which are typically scarce in 
clinical settings. Here, we present a machine learning emulator of the simulation 
model based on a large population of virtual cases, to accelerate re-entry risk pre- 
dictions and circumvent the need for HPC in the clinic. The overall model pipeline 
is presented in Figure 3.6. 


Graph Graph Graph 
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layer layer layer 


Linear} Activation 
layer function 
0 (no re-en 
+ + + H >) ob} “) 
ReLU ReLU ReLU / | Retu 1 (re-entry) 


— ZA 


Cardiac anatomy Graph-based input based on Graph convolutional neural network Classification module Final prediction 
and scar properties the 17 AHA segments 


Fig. 3.6: Proposed pipeline for the arrhythmic risk prediction model using a graph- 
based representation of the heart, including global and local properties. The model is 
based on a graph convolutional network, GCN, and is trained on ground truth labels 
(reentry or no reentry) obtained from numerical simulations. 
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3.2.4.1 Training Data Generation 


Our predictive model is developed using data obtained from the original dataset de- 
scribed in [14]. This dataset comprises 129 cases; 30 original image-based patient- 
specific cardiac models, plus an additional 99 augmented models with varying scar 
size and transmurality, amounting to a total of 1 to 5 cases per patient. Cases origi- 
nating from the same patient are used for either training or testing, not both, to avoid 
training data information leaking into the test set. Given its small size, the dataset is 
split into training and testing sets with a 80%-20% ratio, yielding 104 training cases 
(25 patients) and 25 test cases (5 patients). 

The ground truth label for each case is derived from the results of the numerical 
simulations determining ventricular activation times, reentry circuits, and subse- 
quent risk of arrhythmia in [14]. The binary label extracted from these simulations 
represents the presence (1) or absence (0) of a reentry observed in each simulated 
case. 


3.2.4.2 Data Representation 


Given the complexity of the anatomical models, dimensionality reduction and feature 
extraction are a necessary pre-processing step to capture the data efficiently before 
it can be used as input to the machine learning model. We encapsulate structural 
properties of the cardiac models at a global scale by building a graph-based repre- 
sentation of each model based on the 17 AHA left ventricular segments [23]. Each 
graph consists of 17 nodes and 72 edges; this structure is identical for all cases. Each 
node corresponds to one AHA segment, while edges represent an existing connection 
between neighbouring segments. This representation allows us to implicitly encode 
the structure of the heart in an anatomically accurate manner, instead of having to 
represent this information in tabular format. In addition, we include two regional 
tissue properties computed for each segment: segment volume (ml) and scar volume 
(ml). For segments where no scar is present, the latter value is set to zero. These 
local properties, which vary from case to case, are encoded within a feature vector 
at each node corresponding to the relevant segment. 


3.2.4.3 Model Architecture 


Graph convolutional networks extend on traditional convolutional neural network ar- 
chitectures, originally developed for grid-like Euclidian spaces, by applying the same 
filtering principles to graph-like data [26]. We develop a model using GCN layers 
to process the data encapsulated by our graphical inputs, followed by a classifica- 
tion module to perform binary classification. We propose a GCN-based graph-level 
prediction architecture with the following structure: 


1. One input GCN layer size (2,h) followed by rectified linear unit (ReLU) activation. 
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2. Two hidden GCN layers size (h,h) to convolve the node-wise features, followed 
by ReLU activation. 

3. One output GCN layer size (h,1) to obtain a single node-wise feature. 

4. One fully connected linear layer size (17,1) to obtain a single output feature from 
the 17 node-wise outputs. 

5. A final sigmoid activation layer to obtain the final prediction as a probability 
between 0 and 1. 


Typically, the size h of the hidden layers is chosen as an intermediate between the 
input and output sizes. Given that each node has two input features and one output, we 
select h=2 as an adequate hidden layer size. Only two hidden layers are implemented 
to avoid overfitting, given the small dataset size. The model has a total number of 39 
parameters, less than half the size of the dataset, which also helps to minimise the 
chance of overfitting. 


3.2.4.4 Training and Evaluation 


Binary cross-entropy is chosen as a suitable loss function to train the model, as 
it measures the difference between predicted probabilities and ground truth binary 
labels. Given the small size of the dataset, we train our model for 200 epochs with a 
batch size of 1, a learning rate of 0.001, and the Adam optimiser [27]. 

The model is evaluated using different metrics to reflect its performance in pre- 
dicting presence or absence of re-entry. Accuracy measures the correctly classified 
observations across both outcomes, while F1 score combines precision and recall to 
measure the performance of the model in classifying a single type of outcome. F1 
score is the typical metric used to evaluate the performance of a binary classifier. 


3.2.4.5 Implementation and Software 
We used Python v3.10, with Pytorch v1.4 and Pytorch Geometric v2.3.1. The model 


was trained using Google Colab with an NVIDIA T4 GPU, 12.7GB of system RAM, 
and 78.2GB of disk memory. 


3.3 Results 
3.3.1 Mesh Pipeline 


Figure 3.7 shows a subset of generated meshes of the pipeline. Using the entire 
height of a surface mesh (e.g., Figure 3.2) produces CARP files which are roughly 
twice as large in size compared to those of [14]. 
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Fig. 3.7: 25 out of 30 generated meshes from our meshing pipeline using 50 % 
FWHM. Each mesh contains heart tissue (green) and scar tissue (red). 


In order to reduce the mesh size to similar levels as [14] we decided to restrict 
the mesh generation to the lower 25% part of the biventricular mesh geometry. 
Furthermore, this step also greatly reduces computation times for EP simulations. In 
cases where no scar tissue is present within the 25% lower part of the mesh geometry 
we replace scars using the scar injection algorithm from Section 3.2.2.2. 


3.3.2 EP Simulations 


Despite the reduction of mesh sizes, we were unfortunately not able to complete 
all EP simulations for the entire biventricular population in the allotted time frame. 
Because of this circumstance, we decided to use the data from [14] for the arrhythmic 
risk prediction model. 


3.3.3 Performance of Arrhythmic Risk Prediction Model 


The simulation results taken from [14] are summarised in Figure 3.8. Of the 129 
simulations, re-entry was induced in 53 cases, while the remaining 76 cases resulted 
in a healthy ventricular activation pattern with no re-entries observed. Here, we 
consider cases with induced re-entry as the positive class. After splitting the dataset 
into training and testing sets, we obtain 40 positive cases in the training set (64 
negative cases) and 10 positive cases in the testing set (15 negative cases). We 
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Observe a trend between increasing scar size and occurrence of re-entry, and largest 
percentage of re-entries in areas of the left ventricle near the left descending artery 
(LAD), as seen in Figure 3.8. 
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Fig. 3.8: (left) Distribution of re-entry occurrences in 129 simulated ventricular 
activation sequences. Presence and absence of re-entry are indicated by 1 and 0, 
respectively. The scar location is defined using the AHA segment corresponding 
to the highest scar volume present. (right) Summary of re-entry occurrences and 
frequency in related arteries. 


Our model was trained to perform binary classification on the labelled cases 
provided. The training time for 104 cases was 1 minute and 37s. When evaluated on 
the 25 test cases, the trained model had an accuracy of 88.0% and the F1 score was 
48.0%. The trained model predicted the outcome of a single case in 0.01 seconds. 


3.4 Discussion 


In this project, we have presented a novel approach for biventricular mesh gener- 
ation and rapid arrhythmic risk assessment in post-infarction patients. Our work 
contributes as an improvement of the [14] work in a variety of ways. 

A comprehensive workflow was employed for the processing and preparation of 
cardiac imaging data, with the ultimate goal of creating patient-specific 3D models 
of the ventricles. Several key steps and methodologies were applied to achieve this. 
Automated contour generation for the LV with the software ADAS demonstrated 
high accuracy with minimal manual adjustments, while manual segmentation was 
applied to the RV. Refinement of extracted contours, unified surface mesh generation, 
and the incorporation of realistic cardiac muscle fiber information contribute to the 
creation of a detailed 3D model. 

Various augmentation techniques were introduced to enhance the characteriza- 
tion of scar tissue within the LV. These techniques not only improve the precision of 
scar modeling but also enable the exploration of diverse scar geometries for in depth 
analysis. Firstly, the scar tissue within the LV was characterized using the FWHM 
threshold method. To diversify the scar geometries, three distinct scar configura- 
tions were generated for each LGE MRI scan. Secondly, to further enrich the scar 
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population within each mesh, an innovative artificial scar generation algorithm was 
developed. 

As a result, our innovative mesh generation pipeline successfully produced a 
diverse array of patient-specific biventricular meshes featuring realistic scar geome- 
tries, thereby providing valuable assets for assessing risk of arrhythmia through 
mechanistic EP simulations. Due to time constraints, we weren't able to complete 
the simulations for all biventricular populations. Instead, we used the already present 
data of [14] for the arrhythmic risk prediction model. 

Our predictive model was able to classify reentry events with a high accuracy, but 
a low F1 score, which may be due to the small size of the dataset. The architecture 
of the model was designed to account for this by including relatively few trainable 
parameters. However, the model may still struggle to learn these parameters correctly 
given the limited number of training samples. The small number of input features 
extracted from the dataset may be another potential limitation. We believe that 
the model performance would improve when applied to a larger, richer dataset. 
Additionally, the structure of the proposed model makes it readily extendable to 
predict the location of reentries by performing classification at node-level instead of 
at graph-level. 

Despite its limitations, our model was able to predict the outcome-of-interest 
significantly faster than the original numerical simulations. Typically, one EP simu- 
lation can take upwards of 5 hours to complete for a single patient, using approaches 
and computational resources similar to the ones used in [14]. In contrast, our GCN 
model was able to predict presence of arrhythmia in just 0.01 seconds. This im- 
provement in prediction time represents a speed-up of over five orders of magnitude. 
This is particularly significant in clinical settings, where decision-making may be 
time sensitive and powerful HPC resources may not be available. Furthermore, rapid 
prediction times may also be advantageous in further computational studies, by 
generating large populations of virtual cases with the aim to investigate parameter 
inference and explore trends in cohorts of healthy and diseased patients [28, 29]. 

In conclusion, we presented a novel cardiac meshing pipeline and ML-based 
predictive model to determine risk of arrhythmia in post-MI patients. Our pipeline 
allowed for reconstruction of personalised and physiologically accurate 3D meshes 
of the human ventricles based on 2D LGE MRI scans. Using the patient-specific 
biventricular meshes, we simulated electrical propagation on the ventricles, and 
determined risk of arrhythmia due to reentry. Our predictive model was able to learn 
from the mechanistic simulations and significantly accelerate prediction times. This 
work paves the way towards more accurate and faster computational methods for 
arrhythmic risk stratification, helping to improve therapeutic strategies and long- 
term outcomes of post-MI patients. 
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Abstract During pregnancy, fetal distress requiring clinical intervention can be 
difficult to accurately monitor and diagnose, necessitating technological improve- 
ments to bring clear information to patients and their care teams. Non-invasive fetal 
Electrocardiography (NI-fECG) monitoring may allow for earlier and more reliable 
detection of global cardiac ischemia due to hypoxia. However, the low signal-to-noise 
ratio of the fetal heartbeat relative to the maternal heartbeat remains a challenge. 
To enable reliable recognition of ischemia in NI-fECG, we propose an approach 
that combines simulating a pregnant torso and an unsupervised machine-learning 
method. Three stages of fetal cardiac ischemia: none (healthy), moderate, and se- 
vere, were introduced to the model. For each case, Electrocardiograms (ECGs) were 
simulated with the standard 12 leads, plus 3 additional abdominal leads. Unsuper- 
vised Multiple-Kernel Learning (MKL) with k-means clustering identified changes 
consistent with fetal cardiac ischemia despite noise from the parental heart. Thus, in 
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this early proof-of-concept investigation, our results suggest that NI-FECG may offer 
a means for detecting global cardiac ischemia. 


4.1 Introduction 


The most common techniques for monitoring fetal distress in clinical practice are 
single hand-held Doppler ultrasound or Cardiotocography (CTG) which provide 
information about the Fetal Heart Rate (FHR). However, Doppler ultrasound tech- 
niques cannot provide estimates of Fetal Heart Rate Variability (FHRV) and the 
accuracy of FHR estimates is technically limited by factors such as low sampling 
frequency and imprecise alignment with the fetal heart. Due to the lack of robustness 
of FHR monitoring, Caesarean deliveries are often performed with limited evidence 
of reducing adverse long-term neonatal outcomes [1]. 

Recent studies have highlighted the potential of fetal Electrocardiography (fECG) 
in providing a more reliable assessment of fetal distress. Currently, fECG is obtained 
viaa scalp electrode which is directly attached to the fetal scalp [2]. As this procedure 
is invasive and requires the amniotic membrane to be ruptured, it is not routinely 
performed. Non-Invasive fetal ECG (NI-fECG) provides an alternative approach for 
continuous in-utero fetal monitoring with negligible risk to the parent and fetus, as 
it only requires skin-contact electrodes on the parent's skin to obtain information 
about fetal cardiac activity [3]. Analysis of ECG waveforms is routinely performed 
in clinical practice to provide information about cardiac function in adults. However, 
interpretation of NI-fECG is more challenging due a lack of standardisation of NI- 
fECG acquisition protocols and low signal-to-noise ratios in abdominal ECGs due 
to signal interference from the parental heart, small fetal myocardial volumes and 
parental and fetal motion artefacts [4]. 

Consequently, sophisticated post-processing techniques are required to extract 
the NI-fECG from the abdominal signal mixture. Common NI-fECG extraction 
techniques include adaptive filtering [5, 6, 7, 8, 9], Kalman filtering [10, 11, 12] 
and wavelet transform [13, 14]. Furthermore, blind source separation is a technique 
that leverages principle component analysis (PCA) and independent component 
analysis (ICA) to separate the mixture of independent signals in the abdominal ECG 
[15, 16, 17]. Further NI-fECG extraction techniques include template subtraction 
which involves subtracting a synthetic parental ECG which is generated by using 
parental QRS complex estimates from the abdominal signal to create a template 
parental ECG [18, 19, 20, 21, 22, 23]. While effective for single-channel ECG 
denoising, NI-fECG extraction techniques are significantly limited by the need for 
precise parental QRS complex detection, which is challenging in the case where the 
parental and fetal R wave temporally overlap [24]. 

Deep neural networks have recently emerged as an alternative NI-fECG extraction 
approach [25] with several studies proposing deep convolution encoder-decoder 
network implementations to extract the NI-fECG from the abdominal signal mixture 
[26, 27, 28, 29]. However, there are many challenges of such deep learning based 
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methods. Such challenges include training the model to directly extract the NI- 
fECG due to it being significantly weaker than the parental ECG as well as the 
contamination of the extracted NI-fECG with the residual parental ECG which 
makes identification of the P and T waves challenging. Moreover, the lack of gold- 
standard training data and NI-fECG databases that can be used to evaluate NI-fECG 
extraction algorithms may hinder the robustness and generalisation performance of 
the model[30]. 

An additional approach to analysing ECG signals through the use of artificial in- 
telligence is unsupervised machine learning, which enables the analysis of the natural 
distribution of the data. Due to the complexity of data sets, unsupervised machine 
learning is often used in conjunction with dimensionality reduction techniques to 
allow the number of dimensions of a set of inputs to be reduced while preserving the 
information provided. Principal component analysis (PCA) is the most commonly 
used dimensionality reduction technique [31]. However, when working with clini- 
cal data it can fall short, as it is monovariate and linear. Other techniques, such as 
Multiple Kernel Learning (MKL) have been used to analyse biological signals as a 
multivariate and non-linear alternative [32, 33, 34]. 

In clinical practice, monitoring of the ST segment in adult 12 lead ECGs is com- 
monly performed to diagnose myocardial ischemia [35]. ST segment analysis has 
additionally been identified as a marker of fetal distress due to hypoxia [36, 37]. 
However, use of fetal ST analysis remains limited as it is currently carried out us- 
ing invasive fetal scalp electrodes. Notably, Clifford et al. compared the FHR and 
ST deviation derived from fetal scalp electrode data and non-invasive abdominal 
electrodes, highlighting their clinical similarity [38]. This provides an opportunity 
to leverage computational cardiac electrophysiology to acquire synthetic NI-FECGs 
and develop machine learning techniques to enhance detection of ST segment de- 
viation indicative of fetal ischemia, which could improve clinical decision-making 
regarding when a Caesarean section should be performed to optimise parental and 
fetal outcomes. Herein, we performed fetal and parental cardiac electrophysiology 
simulations using an image-based finite element model of a pregnant torso to acquire 
realistic simulated NI-fECGs in healthy and diseased conditions to determine the 
extent to which we can detect fetal ischemia. 


4.2 Methods 


To generate a data set, 10 simulations were conducted with openCARP [39] followed 
by a machine-learning based approach to characterise ischemia in the NI-FECG. 
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4.2.1 Geometrical Mesh Construction 


In this study, the geometrical mesh of a pregnant torso developed by Uv et al. [4] was 
used. Briefly, a CT-based cardiac biventricular geometry [40] was combined with an 
MRI-based whole body pregnant model [41] to represent the parental and fetal hearts 
embedded in an anatomical pregnant torso. The parental heart was scaled, translated 
and rotated to obtain the fetal heart. The epicardial, myocardial and endocardial layers 
were defined as 30%, 25% and 45% of the cardiac wall thickness respectively [42]. 
Myocardial fibre orientations were assigned to the model using a previously validated 
rule-based method [43]. 


4.2.2 Electrophysiological Modelling 


4.2.2.1 Ionic Models 


The ten Tusscher model of the human ventricular myocyte was used to determine the 
ionic current properties of the parental heart [44]. For the fetal heart, the ionic model 
was modified to reflect increasing ischemic levels in three stages: healthy, moderate 
and severe ischemia. Table 4.1 presents how the default properties were changed 
following the work presented in [45]. For the fetal epicardial, mid-myocardial and 
endocardial layers, the membrane transient outward conductance gro was set to 0.073 
Sm! [3]. 

In the myocardium, the intracellular and extracellular conductivities were set to 
gi = (0.27, 0.081,0.045) Sm! and ge = (0.9828, 0.3654, 0.3654) Sm! in the 
longitudinal, transverse and normal directions respectively and the torso was given 
an isotropic conductivity of 0.216 Sm7! [3]. Tissue propagation was modelled with 
the monodomain approach. 


Table 4.1: Summary of the parameters of the ten Tusscher ionic model of the human 
ventricular myocyte [44] modified to simulate healthy, moderate and severe ischemia 
in the fetal endocardial layer. 


Cell Model Parameters Healthy Moderate Ischemia Severe Ischemia 
[K* ]o (mmol/l) 54 12.5 15.0 

ENa» £ca,L (9) 100 75 50 

Pnak (%) 100 100 30 

knaca(%) 100 100 20 

Vmaxup (%) 100 100 90 


V, o4 (1) 100 100 5 
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4.2.2.2. Stimulation Protocol 


Both the fetal and parental heart were stimulated following the protocol proposed by 
Durrer et al. [46]. For all simulations, the parental heartbeat was set to a cycle length 
of 667 ms, which corresponds to a healthy adult heartbeat during the second trimester 
of 90 bpm [47]. For the fetal heart, simulations for the following different heart cycle 
lengths (heart rates) were conducted: 462 ms (130 bpm), 400 ms (150 bpm), 353 ms 
(170 bpm). For a gestational age over 32 weeks, these heart rates correspond to the 
range of normal fetal heart rates between 110—160 bpm [48]. 


4.2.3 Extracellular Potential Measurements 


The extracellular potentials were recovered in a post-processing step with openCARP 
at the electrode locations defined in Figure 4.1. The four limb and six precordial 
electrodes match the standard electrode positions to calculate the 12 lead ECG. 
Figure 4.1 also illustrates the calculation of three abdominal leads derived from 
electrode positions from previous work [4]. 

Since the extraction of the NI-fECG without distorting its waveform morphology 
is an unsolved problem, an ECG showing only fetal activity was required. Hence, we 
assumed that the simulation results with mixed activity from the parental and fetal 
heart could be obtained approximately by superimposing simulations results where 
only one heart was active at a time. To obtain a pure fetal signal, simulations where 
only the parental heart was active were performed and subtracted from simulation 
results with parental and fetal cardiac activity. 

To investigate how sensitive our proposed analysis protocol was to inadequately 
extract NI-fECGs, we developed a method to linearly suppress the parental activities 
in the abdominal recordings. Accordingly, chest and abdominal leads were scaled 
and subtracted from each other, using only leads that captured the parental activity 
in the same direction, e.g. leads A52 and aVF from the 12 lead ECG. 


4.2.3.1 Simulation Setup 


By changing the fetal heart rate and inducing three stages of increasing ischemia 
in the fetal heart, nine simulation sets of 5s recordings were generated. To reduce 
overall computation time, a simulation where only the parental heart was active was 
conducted to retrieve the pure fetal signal for all nine cases. 
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Fig. 4.1: Visualisation of the pregnant torso model illustrating the electrode positions 
and the calculated abdominal leads. 


4.2.4 ECG Signal Analysis 


Since NI-fECGs are complex signals with multiple sources (leads) aligned in time, we 
applied unsupervised machine learning analysis to a dimensionality reduced space 
to determine whether we could differentiate between healthy and ischemic cases in 
the NI-fECGs. Using an analysis pipeline adapted from the work by Jimenez-Perez 
et al. [49], we firstly used a delineation deep learning network [50] to identify the 
QRS-complexes, T-waves, and S T-segments in all the leads of the different simulated 
NI-fECGs. Secondly, we applied MKL on the data set to generate the dimensionality 
reduced space. Finally, we perform a clustering analysis, characterising every cluster 
to determine the discriminative power of the method. 


4.2.4.1 Signal Delineation 


The used ECG delineation network enabled the segmentation of biomarkers that 
we considered relevant for this problem. Separating the NI-fECG from the parental 
ECG was outside of the scope of this project. However, since our goal was to analyse 
the NI-fECG morphology, we applied the delineation network on the pure NI-fECG 
which only shows fetal cardiac activity. Following this, the mask was applied to the 
original signal in which both fetal and maternal hearts were active. Therefore, the 
ideal segmentation enabled the focus of the analysis to remain on the discrimination 
capacity of our solution. 
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4.2.4.2 Dimensionality Reduction 


To increase the number of cases in the dataset, we treated each cardiac cycle as a 
separate record. Since the simulations used biventricular geometries, no P wave was 
generated. Hence, the start of the cardiac cycle was taken as the beginning of the 
QRS complex and the end of the cardiac cycle was marked at the end of the T wave. 
All segments were then resampled to be of the same length and aligned to allow the 
MKL algorithm to generate the reduced dimensionality space. 


4.2.4.3 Cluster Analysis 


K-means clustering was used to partition the data into clusters to determine if under- 
lying differences between the healthy, moderately ischemic and severely ischemic 
NI-fECG could be detected. The number of clusters was determined using the sil- 
houette method [51], which allowed us to choose the number of clusters that best 
separates the data. 


4.3 Results 
4.3.1 Simulation Results 


From each simulation a 15 lead ECG was constructed containing data from the 
standard 12 lead ECG with 3 additional abdominal leads: A41, A52, and A63 
(Figure 4.2). Ideal NI-fECGs were constructed by subtracting the known parental 
contribution from each lead. However, the actual raw signals from our model include 
substantially more contribution from the parental heart than the fetal heart, obscuring 
even dramatic changes to the ST segment and T wave. Application of a simple linear 
method to suppress the signal from the parental heart enabled the visualisation of 
ischemic changes in the NI-fECG signal, although this still failed to reproduce the 
ideal signals (Figure 4.3). 


4.3.2 Signal Analysis 


The delineation network was applied on the pure NI-fECG signal and the mask of the 
onsets was applied on the original signal, which is contaminated with the parental 
ECG. Results for a single lead are shown in Figure 4.4. The biomarkers chosen for 
the analysis were the QRS complexes (in red), T waves (in blue) and the ST segments 
(the interval between the end of each QRS complex and the beginning of the T wave). 

Each cardiac cycle was considered as a single input for the MKL method. After 
aligning the cycles, each segment was interpolated to a fixed value of samples, 
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Fig. 4.2: Simulated 15 lead ECG of a healthy fetus with a heart rate of 130 bpm. 
Note, the fetal contribution is by the parental heartbeat, even in the abdominal leads 
A41, A52, and A63. Arrows indicate the beginning of each fetal heartbeat. 
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Fig. 4.3: Example of a healthy, moderately ischemic, and severely ischemic fetal 
heartbeat at 130 bpm. The parental heartbeat was removed ideally (left) and by a 
simple linear suppression technique (right). Arrows indicate the beginning of each 
fetal heartbeat. 


allowing the database to be stored in an array. Finally, the MKL algorithm is applied 
on the data set to generate the output space. This process was applied to the data 
set of ideal abdominal ECGs where only the fetal heart was active as well as to 
the original abdominal ECG where the parental cardiac activity was suppressed. 
Figure 4.5 shows the output space for each case following dimensionality reduction. 

The output of the MKL algorithm has N — | dimensions, where N represents 
the number of inputs or records. However, the highest variability is concentrated in 
the first dimension. For the clustering analysis, we considered the dimensions that 
accumulated 95% of the variability of the output space. 

To cluster these results, the output space served as an input for the k-means 
algorithm. The number of clusters was decided by using the silhouette method 
leading to an optimum number of clusters between 3 and 16. From the simulations, 
it was known whether a data point corresponded to healthy, moderately ischemic, or 
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Fig. 4.4: Segmented ECGs for the ideally isolated NI-fECG (left) and the original 
NI-fECG (right) which includes the parental and fetal cardiac activity. The segments 
for the QRS complex (red) and the ST segment (blue) are detected in the isolated 
NI-fECG (left) and are applied to the original NI-fECG (right), both from the aVF 
lead in this case. 
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Fig. 4.5: Output spaces for the different ischemic levels from the abdominal ECG 
leads showing only fetal cardiac activity (left) and showing the original ECG after 
applying the parental suppression method (right). The markers indicate the cluster 
to which each point belongs, revealing four different clusters detected for the data 
set displayed on the left and three clusters detected for the data set on the right. 


severely ischemic data. Using this information, Figure 4.5 illustrates the distribution 
of the different cases in the output space, highlighting the ability of the proposed 
method in distinguishing the three levels of ischemia. For each level of ischemia, the 
total number of records within each cluster was also computed, as summarised in 
Figure 4.6. From this, it is evident that healthy data points are mostly assigned to the 
same cluster. Data based on the original abdominal ECG with parental suppression 
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Fig. 4.6: Summary of the distribution of the data points across all clusters based 
on the abdominal lead recordings with only fetal cardiac activity (Left) and on the 
original abdominal ECG after applying the parental suppression method (Right). 
Each coloured bar indicates the percentage of data points assigned to a specific 
cluster for each ischemic level. 


is more likely to be assigned to the cluster that is mostly associated with healthy data 
points than data based on the pure NI-fECG. Differentiating between moderately 
and severely ischemic NI-fECGs based on the cluster allocation is ambiguous for 
both data sets. 


4.4 Discussion 


The proposed model framework successfully generated synthetic abdominal ECGs 
containing signal contributions from both the parent and fetus. Moreover, a simple 
linear NI-FECG extraction method was implemented to suppress the parental signal, 
allowing dramatic ST segment changes to be identified in the ischemic fetus. In 
addition, the feasibility of using common machine learning techniques such as MKL 
for dimensionality reduction and k-means clustering to analyse the ECG signals was 
demonstrated, providing qualitative and quantitative evidence of distinct clusters 
describing the degree of fetal ischemia. 

Due to the high computational cost of performing whole torso electrophysiol- 
ogy simulations, only monodomain simulations were performed in this preliminary 
investigation. Furthermore, the torso was modelled as a uniform medium hence vary- 
ing tissue conductivities that affect the wavefront propagation were not considered. 
Performing bidomain simulations and accounting for additional tissue types such as 
bone, lungs and the fetal vernix caseosa would provide more realistic conduction 
throughout the torso and thus, more realistic ECGs. Notably, the Paniflov ten Tuss- 
cher ionic model used in this study was adapted to model fetal ischemia using data 
based on an adult human heart. Modifying the ionic model to reflect the ischemic 
fetal heart would be important to obtain more clinically realistic NI-fECG signals. 
Moreover, as this was a preliminary investigation, we focused on modelling severe is- 
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chemia only in the endocardial layer of the fetal heart to ensure the ischemic changes 
were sufficiently large to detect a change in the NI-fECG. Future studies simulating 
lower levels of ischemia in the whole fetal heart would present a more clinically rele- 
vant scenario in which the fetus would have a greater chance of resuscitation. While 
the pipeline developed in this study presents a promising proof-of-concept, using a 
simple linear method to suppress the parental heart signal is not realistic for clinical 
settings. NI-FECG extraction remains a challenging problem, therefore as NI-fECG 
extraction algorithms continue to evolve and improve in the future, techniques should 
be validated using real clinical data and incorporated into the workflow developed 
in this study to enhance the robustness of this proposed method. 


4.5 Conclusions 


We developed, to our knowledge, a novel workflow to detect fetal ischemia by 
simulating parental and fetal cardiac electrophysiology in a pregnant torso model 
and using MKL with k-means clustering to identify ischemic changes in signals with 
amplified fetal contribution. With continued development and further optimisation, 
our methodology may enable improved detection of fetal cardiac ischemia to better 
inform clinical decision-making. In future work, the pipeline will be applied to 
additional pregnant torso models with lower levels of ischemia applied globally 
to the fetal heart to better characterise fetal distress non-invasively. Additionally, 
the impact of incorporating time-dependent factors such as fetal movement and 
physiological artefacts due to parental respiration and muscle contraction should be 
explored. 
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Reconstruction of a Pancreatic Beta Cell 
Network From Heterogeneous Functional 
Measurements 


Roshni Shetty, Radhika Singh-Agarwal, Stefan Meier, Christian Goetz, Andrew G 
Edwards 


Abstract Intercellular heterogeneity is fundamental to most biological tissues. For 
some cell types, heterogeneity is thought to be responsible for distinct cellular pheno- 
types and functional roles. In the pancreatic islet, subsets of phenotypically distinct 
beta cells (hub and leader cells) are thought to coordinate electrical activity of the 
beta cell network. This hypothesis has been addressed by experimental and compu- 
tational approaches, but none have attempted to reconstruct functional specialization 
directly from measured heterogeneity. To evaluate if electrophysiologic heterogene- 
ity alone can explain these specialized functional roles, we created a population 
of human beta cell models (via genetic algorithm optimization) recapitulating the 
heterogeneity in an extensive patch clamp dataset (1021 pancreatic cells). Then we 
applied the simplified Kirchhoff network (SKNM) formalism to simulate activity 
of that population in a connected beta cell network. We could not immediately ob- 
serve cells with obvious hub or leader phenotypes within the heterogeneous network. 
However, with this study we built the basis for further ’ ground-up” investigation of 
relationships between beta cell heterogeneity and human islet function. Moreover, 
our workflow may be translated to other tissues where large electrophysiologic data 
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sets become available, and heterogeneity is thought to influence tissue function e.g. 
human atria. 


5.1 Introduction 


Intercellular heterogeneity is the cornerstone of functional complexity and enables 
coordinated responses and specialized functions of cells within a tissue [1]. Hence, 
it has become increasingly important to study intercellular heterogeneity to unravel 
the mechanisms underlying the basic function of many organs, as well as disease and 
development. However, the ability to make detailed assessment of cellular properties 
in a sufficient number of cells to properly characterize heterogeneity is often limited 
by experimental throughput. Historically this has certainly been true, and remains 
true, for electrophysiological characterization of many tissues (the brain being an 
important exception). It is only rarely that characteristics other than the Gaussian 
mean and some measure of distribution spread are reported for most electrophysi- 
ologic variables. Expression measurements at the RNA and protein levels are often 
provided in more detail, such as for markers measurable by flow cytometry. However, 
this is also rare for ion transporters, the building blocks of cellular electrical activ- 
ity. To date, these challenges in experimental data collection have largely prevented 
modeling approaches that have tried to systematically link heterogeneity in cellular 
properties to macroscopic electrophysiologic function. 

The pancreas is an organ with both endocrine and exocrine functions. The en- 
docrine cells (alpha, beta, gamma, delta, and epsilon) are arranged in small clusters 
(several hundred to one thousand cells) called islets, which serve as mini-organs. 
Cells of an islet work together to maintain blood glucose homeostasis [2]. Beta cells 
that synthesize and secrete insulin, account for approximately 70% of human islet 
cells [3]. Experimental studies of isolated beta cells show that they exhibit high 
intercellular variability in terms of the insulin the secrete in response to both basal 
glucose concentration (2-3 mM), and upon increasing glucose to simulate prandial 
conditions. This functional heterogeneity is of great importance to normal islet be- 
havior and regulation of insulin release [4]. A prominent hypothesis in the field of 
beta cell biology is that subsets of beta cells served different roles in determining 
both electrical activity and insulin release (e.g. leader, hub, and first response cells) 
[5, 6]. However, the metabolic and electrical profile of these different integrated phe- 
notypes is yet to be definitively characterized. Here, we attempt to answer whether 
experimental data describing heterogeneity in a range of cellular electrophysiologic 
properties can explain the specialized functional roles of heterogeneous beta cells 
within an islet. 

In this study, we combine in vitro data obtained from a recent and comprehen- 
sive patch clamp data set of isolated human beta cells [7] with a computationally 
designed framework for reconstructing that heterogeneity in a population of in silico 
beta cell models (PoM) approach [8, 9]. In this approach, key model parameters (e.g., 
channel conductances, transporters maximal transport velocity) are varied to reca- 
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pitulate measured reference data. To align the simulated population distribution with 
the experimental distribution, we employ genetic algorithms as an able optimization 
method to find appropriate parameter adjustments. These are "good enough" solu- 
tions to the optimization problem, but are generally not unique. Subsequently, the 
diverse simulated population is coupled in a tissue-level architecture, with heuristics 
to specify intercellular coupling in a manner that is consistent with measurements, 
but random with respect to the cellular electrical phenotype. This in silico islet can 
then be perturbed to gain insights into the functional significance of heterogeneity 
at a tissue level. 


5.2 Methods 
5.2.1 Data and Model 


The experimental data utilized in this study was obtained from Camunas-Soler et al. 
[7] (https: //github.com/jcamunas/patchseq) which comprises an extensive 
collection of patch clamp data encompassing 1021 pancreatic cells. Within this 
data set, information concerning six key electrophysiological outcome metrics were 
extracted from a series of 4 distinct voltage clamp protocols. The metrics were: (1) 
early (EE) and (2) total exocytosis (TE) (in fF/pF); (3) early and (4)late Ca?* current 
(in pA/pF); (5) peak Na* current (in pA/pF); and, (6) Na* channel half-inactivation 
(Vna; r,in mV). These data were measured in 180 healthy beta cells. Figure 5.3 shows 
the experimental distributions of the outcome metrics. To implement the voltage 
clamp protocols and derive the six outcome metrics for each member of the in silico 
population of models, we used MATLAB (release 2022B). These simulations were 
executed on the human pancreatic beta cell model developed by Riz et al. [10]. 


5.2.2. Workflow 


Figure 5.1 shows the workflow that was followed to match the simulated model data 
distributions to the experimental electrophysiological data distributions [7] which 
is computationally expensive. Subsequently, the optimized parameters were used to 
create a 3D heterogeneous model of a human pancreatic islet network. First, we em- 
ployed a sensitivity analysis to identify the most important parameters related to Na* 
-and Ca?* outcomes (i.e., EE and TE (in fF/pF), early and late Ca?* current peaks 
(in pA/pF), peak Na* current (in pA/pF), and Vhalf (in mV) to reduce the parameter 
space and increase computational efficiency. From the initial 86 parameters, only 
eight were selected for optimization. Subsequently, an iterative process of model 
generation and parameter optimization through a genetic algorithm was employed to 
map the model distributions to the experimental distributions. Once the fits matched 
sufficiently, several populations were generated. Finally, these populations were em- 
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bedded in a network of connected beta cells forming a heterogeneous 3D model of 
a pancreatic islet. 


Fig. 5.1: Workflow for parameter selection, optimization, model generation and 3D 
simulations. 


5.2.3 Sensitivity Analysis 


A sensitivity analysis was performed to better understand the impact of individual 
parameters on the Na* -and Ca?*-related outcome parameters, and thereby constrain 
the search space for the optimization algorithm. A simple linear scaling —by factors 
of 4, 2, 1, 0.5, and 0.25 — was applied to each model parameter individually, while 
the remaining parameters were held at their baseline value. First, the voltage-clamp 
protocols were pre-paced for at least 10s after which the protocols were evaluated. 
The changes in outcome parameters for each scaling factor were compared to their 
baseline value. 


5.2.4 Genetic Algorithm Optimization and Population Generation 


From the sensitivity analysis, the parameters with most influence on the electrophys- 
iological metrics were identified and collectively defined as set (P). An initial value 
of standard deviation for every parameter belonging to the set (P) was obtained from 
independent small-sample data sources [11], and used as the initial guess of variation 
for these parameters within the beta cell population. For each iteration of the GA, 
a new population of beta cells was generated by applying the standard deviation in 
set (P) to sample a from a log-normal distribution using the baseline value of the 
Riz et al. model as the centre. 750 models were generated in each population. 

To address the bimodal distribution observed in the Na* current half-inactivation, 
we introduced a second Na* current component to the Riz ef al. model. This sec- 
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ondary current had a lower average half-inactivation voltage, resulting in a population 
where each cell expressed either the default high Na* current (as in the original Riz 
et. al model) or the newly introduced low Na* current. The parameter Frachigh gov- 
erned the ratio of models expressing the high and low Na* current components (each 
model cell was assumed to only exhibit one dominant form of the current). An initial 
value of 0.15 was set, resulting in 113 models representing the high Na* current 
and 637 models expressing the low inactivation Na* current. Next, the six outcome 
metrics were assessed for each model in the population. The primary objective of 
the optimization was to closely align the simulated population distributions of these 
metrics with the corresponding experimental distributions. 

Given the non-normal nature of the distributions, parameters were calculated to 
capture skewness (S), kurtosis (K), and the mean-to-standard deviation ratio (u/o) 
for both experimental and simulated distributions of skewed Gaussian metrics (such 
as early and total exocytosis, early and late Ca?* current peaks, and peak Na* current). 
These calculations were facilitated using built-in MATLAB functions, and the GA's 
cost function aimed at reducing the error in (S), (K), and (u/o") between the two 
distributions (Equation 5.1 and 5.2), 


Cost function = arg min E; (k) (5.1) 
E\(k) = 9 ISe- Ss| [Ke — Ks| + |Me/oe — us[os| (5.2) 
k 


where k is the type of metric distribution and subscript e and s denote the value 
derived from the experimental and simulated distribution respectively. 

For the bimodal distribution of the Na* current half-inactivation, a double gaussian 
curve was fitted using the MATLAB curve-fitting toolbox. Parameters u1, 42, 01, 02, 
and the ratio of the two gaussian peaks a1/a» were selected for optimization. Simi- 
larly, the differences between these parameters were minimized in the cost function 
for both experimental and simulated distributions (Equation 5.3 and 5.4). 


Cost function = arg min E(k) (5.3) 


Ex(k) =% la - Hi| + duo, = H| doi, = 01,] + lo, 7 0,1 
k 


|a, /a5, — a1,/a2, | (5.4) 


GA optimization was performed with the Global Optimization Toolbox in MAT- 
LAB [12, 13]. Figure 5.1 outlines the main steps in this process. The standard 
deviations of parameters in set (Pj and the Fracpign were iteratively adjusted in 
each generation to minimize the error between the simulated and experimental dis- 
tributions. The optimized standard deviations for set (P) and FraChign values were 
selected to construct the final population, which was again constructed by log-normal 
sampling of parameters in set{P} as well as Fracyign. For other non-optimized 
parameters (not-identified by the sensitivity analysis), independent measures of in- 
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tercellular variability (standard deviation) were taken from the source data used to 
construct the baseline Riz model, or similar human beta cell datasets ([11]). As 
earlier, the mean of this log-normal variation was always centered at the baseline 
value of the Riz et al. model parameter. Table 5.1 presents a list of all parameters 
varied, along with the corresponding means and standard deviations (normalized to 
the baseline Riz et al. model value). 


5.2.5 3D Beta Cell Network Simulations 


Multiple cell-based models have been developed in recent years [14, 15, 16, 17, 18, 
19], to capture tissue electrodiffussion in varying degrees of detail. A recently pub- 
lished model aiming at balancing the need for cell-level accuracy and computational 
efficiency is the KNM [20]. This model incorporates individual cells with explicit 
representation of the local extracellular space. The fundamentals of this model are 
related to Kirchhoff’s current law which states that the sum of all currents flowing 
into one node (cell) is equal to the sum of all currents flowing out of that node (cell). 
While incorporating individual cells, the computational efficiency of the KNM is 
still comparable to to models in which extracellular and intracellular space are not 
geometrically distinct e.g. the bidomain (BD) model [20]. The KNM can be reduced 
to the simplified KNM (SKNM) by assuming equal anisotropy of intra- and extracel- 
lular conductances, which further improves efficiency with little penalty on realism 
[21]. In this study, we apply the SKNM to study the effect of cellular heterogeneity 
on the integrated function of a connected network of islet beta cells. We assumed 
the extracellular potential to be negligible (equal to 0) which leads to the following 
equations [21] (or pr^ >> Gijj,k): 


dv 1 "P 
Car ~ AK 2, 6/9 VE o tat. (5.5) 
m jeNy 
ds* 
dt Fe(s“, v“) (5.6) 


Cm describes the specific membrane capacitance (in pF/cm?), A* is the mem- 
brane area of cell k (in cm”), Nx is the number of connected cells to cell k, Gi Kis 
the intracellular conductance between cells j and k (in mS), v* is the membrane po- 
tential of cell k (in mV) (analogue for v/) and I 2 is the ion current density through 
ion channels, pumps and exchangers on the membrane of cell k (in pA /cm?). Lastly, 
F, is a function describing the dynamics of a number of state variables s* that model 
the membrane dynamics of cell k [21]. 

For the computational domain we used a human islet with 257 individual beta 
cells reconstructed from confocal microscopy data containing information about the 
positions, radii and contacts between all beta cells, based on ([22, 23]). To define 
which cells are electrotonically coupled via gap junctions, and to set the conduc- 
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tance for those connected beta cells, we computed the overlap (or separating distance) 
between adjacent cells in the confocal dataset. We assumed the gap junction conduc- 
tance to be proportional to that distance, with an average gap junction conductance of 
G mean = 2 X 1077 mS [24]. The individual gap junction conductance Gi K between 
cell j and k is then calculated based on the distance //-* between cells j and k 
normalized to the mean distance lmean between all cells: 


ES 


p = ' Gmean (5.7) 


Imean 

As a baseline, we performed simulations of a homogeneous network i.e., each beta 
cell in the geometry represents an instance of the default Riz et al. model with original 
parameters. To study the impact of electro-metabolic heterogeneity we replace each 
of those default model instances with instances taken from the final population of 
the Riz et al. models. Thus, each cell in the heterogeneous beta cell network has an 
instance of the Riz et al. model with some variation (from default) in the underlying 
parameters (table 5.1), although the distribution across all cells remains fitted to the 
experimental distributions. Finally, the equation system is solved for each time step 
in two individual steps. First, the ordinary differential equations (ODE) of the cell 
model are solved by using the forward Euler method. The time step for solving the 
ODEs is typically ten times higher than for the partial differential equations (PDE) to 
ensure achieving a steady-state in the membrane model. Second, the PDEs are solved 
and the membrane potentials are updated. Parameters for the temporal discretization 
and numerical solver are shown in the supplementary material 


5.3 Results 


5.3.0.1 Sensitivity Analysis 


As anticipated, the parameters related to Ca?* channel and Na* conductance and 
gating primarily affected the outcome metrics of interest (i.e., EE, TE, Vhalf, peak 
Ina, peak Ica, and late Ica), while the parameters related to K* channels did not. The 
following parameters were selected based on the 2- and 0.5-fold scaling results for 
the optimization process by the GA: gna, £caPQ. £CaL. £NaLow; V mNa» VhNa, V mNaLow: 
and VhNaLow (Figure 5.2). 


5.3.1 Genetic Algorithm Optimization and Population Generation 


Two sets of GA optimizations were executed as outlined in Subsection 5.2.4. The 
initial set concentrated on optimizing the standard deviations of gcapg and gcal, while 
minimizing E R,, pertaining to four distributions namely, EE and TE, and early and 
late Ca?* current peaks. The second optimization was for the standard deviations 
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Fig. 5.2: Sensitivity analysis. The top panel shows the results of 2-fold parameter 
scaling and the bottom panel shows the results of 0.5-fold parameter scaling. 


Of £Na. £Nalows WmNA> VhNa» VmNaLow» WhNaLow and the parameter Fracpjgn. In this 
case, ER, was minimized for the peak Na* current, while ER; was minimized for 
the Na* current half-inactivation (mV). 

A detailed compilation of all parameters that were varied in the populations is 
presented in Table 5.1, encompassing the 9 parameters that were optimized us- 
ing GA, as well as the additional 13 parameters whose standard deviations were 
obtained from independent data sources. The ultimate population was generated 
through log-normal variation, using the optimized standard deviations of parameters 
and Fracpign (determined by the GA), as outlined in Subsection 5.2.4. Figure 5.3 
provides comparison between the experimental and simulated metric distributions 
for all outcome metrics. 
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Table 5.1: Parameters varied in the model to create the final population, and their 
corresponding mean (u) and standard deviation (c); (Parameters optimized by the 
GA are in bold font) 


Parameter H c Reference 
VGK max 1.31 1.15 [25 
KGK 1.48 1.44 [26 
VPFKmax 1.00 0.03 [27] 
KpFK 1.00 0.09 [27] 
hprK 1.02 0.14 [27 
VGAPDH,ax 1.02 0.29 [28 
SKv 1.08 0.40 [29] 
&BK 1.26 1.14 [29] 
£caL 1.01 0.15 GA 
2caPQ 1.20 0.68 GA 
&CaT 1.07 0.44 [29] 
SKATP 1.47 1.53 [30] 
SHERG 1.05 0.34 [31 
VnNa 1.01 0.18 GA 
NhNa 1.00 0.12 [29] 
£Na 1.54 1.72 GA 
VinNa 1.00 0.25 GA 
ViNatow 1.02 0.24 GA 
Uh Nato, 0.98 0.16 [7] 
ZNatow 1.33 1.01 GA 
VinNatow 0.97 0.35 GA 
Frachigh 0.21 = GA 


5.3.2 3D Beta Cell Network Simulations 


To assess the impact of heterogeneity on the behavior of an intact beta cell network, 
we compared a network of homogeneous beta cells with a heterogeneous network 
in which cells were sampled from the population optimized to fit the experimental 
distributions. For both simulations, the islet structure we used consisting of 257 beta 
cells with identical connectivity. Figure 5.4 shows membrane potential frames from 
these simulations. In the left column, we see that all beta cells of the homogeneous 
network depolarize simultaneously from a resting membrane potential (RMP) of 
—68 mV to a peak membrane potential of —11 mV. This is to be expected, given 
that all cells in the geometry have identical cellular properties. In contrast, in the 
heterogeneous network no coordinated depolarization can be observed. While some 
beta cells depolarize to a peak membrane potential of 16 mV, most beta cells remain 
at a membrane potentials between —73 mV and —30 mV. 

To explore the role of specific subpopulations in the heterogeneous network, 
we analyzed the time courses of the membrane potentials of individual beta cells. 
Figure 5.5a illustrates the membrane potential of a typical beta cell in the homo- 
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Fig. 5.3: Probability density histograms for the outcome metrics of the experimental 
and simulated population distributions [Blue: experimental distribution (n = 180); 
Red: optimized simulated distribution (n = 750); Purple: overlap] 


geneous network within two periods of bursting behavior (40s - 74s and 185s - 
219s) and a peak membrane potential of —11 mV. In the heterogeneous network 
we identified all beta cells with bursting behavior and evaluated if neighboring cells 
exhibited coordinated activity. Figure 5.5b shows a subset of six neighboring cells in 
the heterogeneous network in which interdependencies could be observed. The most 
significant bursting behavior can be observed for beta cell 141 with peak membrane 
potentials of -27 mV. However, the RMP of all beta cells of the heterogeneous 
network is significantly higher than in the homogeneous network with membrane 
potentials of —50 + 5 mV. To evaluate if beta cell 141 exhibits hub-like behavior, 
we present the behavior of cells that are in its subnetwork. Figure 5.5c shows the 
connections of beta cell 141 with neighboring beta cells. It can be observed that 
the amplitudes of the beta cells diminish with increasing distance from cell 141 - 
see membrane potentials of beta cells 141, 156, 157 and 154 Figure 5.5b. However, 
beta cells 158 and 153 show significantly smaller amplitudes than beta cell 157 
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Fig. 5.4: Membrane potentials at different time steps of a human islet with homoge- 
neous (left column) and heterogeneous cellular properties (right column). 
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Fig. 5.5: Time courses of the membrane potential of individual beta cells in the 
homogeneous (a) and heterogeneous islet network (b) and schematic subset of the 
cell-cell contacts of the islet network (c). 


although they are all connected to beta cell 156. To summarize, only one subset of 
six neighboring beta cells of the heterogeneous network could be observed that show 
coordinated depolarizations. Moreover, the amplitudes of depolarizations decrease 
with increasing distance to the beta cell with the highest amplitudes. Together these 
analyses demonstrate that the higher activity of cell 141 is not sufficient to overcome 
the electrotonic load of its coupled adjacent beta cells. 


5.4 Discussion 


In this study, we have provided a workflow to create a heterogeneous model pop- 
ulation that reflects experimental data to subsequently study the effects of this het- 
erogeneity on a coupled pancreatic beta cell network. The sensitivity analysis was 
the first step to reduce the parameter space and to speed up the model fitting and 
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optimization procedure, while getting an idea about how the model responds to 
changes in parameters. In this study, we used a linear scaling technique similar to 
Meier et al. [32] to scale individual parameters which is rudimentary in comparison 
to Monte Carlo and Bayesian methods [33]. The parameter selection was based on 
the 2-fold and 0.5-fold scaling, because these scalars resulted in a robust change 
without producing unphysiological behavior. Figure 5.2 also shows that the param- 
eters behave in a non-linear fashion. For example, a 0.5-fold decrease in Vnwarow 
caused a more than 2-fold increase in peak Ica and late Ica, while increasing the 
parameter by 2-fold did not meaningfully change peak Ic, and late Ica. Nonetheless, 
the identified parameters obtained by our linear sensitivity analysis provided a good 
starting point for the GA and substantially increased computational efficiency by 
reducing the parameter space to only 9 parameters. 

A major initial assumption was that through implementing a log-normal varia- 
tion within the parameter space would yield skewed distributions for early and late 
exocytosis, early and late Ca?* current peaks, and the peak Na* current, along with 
a double Gaussian distribution for the Na* current half-inactivation. Our findings 
validate this assumption by demonstrating that the optimized final simulated pop- 
ulation closely aligns with the experimental distributions (depicted in Figure 5.3). 
This outcome corroborates the insights from several previous studies [34, 35]. With 
the coherence observed in the simulated distributions, it is credible to claim that the 
simulated population mirrors the inherent heterogeneity observed in experimental 
electrophysiological metrics. 

Through our SKNM network simulations, we showed that randomly applying 
experimentally constrained electrophysiologic heterogeneity tends to stabilize the 
overall behavior of a coupled network of human beta cells. Moreover, we could not 
observe locally specialized functional roles within the network based on intercellular 
heterogeneity. Although there were clear examples of locally important electrotonic 
synchronisation, the electrotonic load was too great to permit broader influence and 
regenerative local or global activation. The properties of neighbouring cells were 
important in determining this local electrotonic load. This can be seen in Figure 5.5b 
and Figure 5.5c, where the electrotonic distance (number of cells from source cell 
141) was less important than the coupled cell's properties in determining their re- 
sponse.Thus, it is likely that either non-random distribution of cellular properties, 
or perhaps reduction or regionalisation of intercellular coupling may be required to 
observe typical hub or leader cell roles. These relationships require further assess- 
ment, and can now be studied in the future by extracting subpopulations from the 
experimental data set and by applying these parameter sets in specific regions of the 
islet. 


5.5 Conclusion 


This study showed a general approach to introduce physiologically realistic hetero- 
geneity in a beta-cell network model by combining a population-of-models approach, 
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with a cell-based tissue electroconduction model - the SKNM. Since randomly ap- 
plying heterogeneity (obtained from the fitted population) to the network model did 
not immediately yield prototypical glucose responses or experimentally observed 
functional specialization within the network, further investigation is required to de- 
termine the potential for regional similarities in cellular phenotype or intercellular 
coupling to permit those behaviors. 

This work defines an initial framework for generating heterogeneous models of the 
human islet that are driven by recent high-volume experimental datasets. Models of 
this type are likely to be useful in understanding the role of intercellular heterogeneity 
in driving specific aspects of network behavior. 


References 


1. Steffen Walherr. Estimation methods for heterogeneous cell population models in systems 
biology. Journal of Royal Society, October 2018. 

2. Kara Rogers. islets of Langerhans, volume 2. Encyclopedia Britannica, 2022. 

3. Joseph W. Kim, John Z.Q. Luo, and Luguang Luo. Chapter 10 - bone marrow mesenchymal 
stem cells as a new therapeutic approach for diabetes mellitus. In Xiao-Dong Chen, editor, 
A Roadmap to Non-Hematopoietic Stem Cell-based Therapeutics, pages 251—273. Academic 
Press, 2019. 

4. Richard KP Benninger and Vira Kravets. The physiological role of B-cell heterogeneity in 
pancreatic islet function. Nature Reviews Endocrinology, 18(1):9—22, 2022. 

5. Gabriela Da Silva Xavier and Guy A. Rutter. Metabolic and functional heterogeneity in 
pancreatic B cells. Journal of Molecular Biology, 432(5):1395—1406, 2020. Islet Biology in 
Type 2 diabetes. 

6. Natalie R. Johnston, Ryan K. Mitchell, Elizabeth Haythorne, Maria Paiva Pessoa, Francesca 
Semplici, Jorge Ferrer, Lorenzo Piemonti, Piero Marchetti, Marco Bugliani, Domenico Bosco, 
Ekaterine Berishvili, Philip Duncanson, Michael Watkinson, Johannes Broichhagen, Dirk 
Trauner, Guy A. Rutter, and David J. Hodson. Beta cell hubs dictate pancreatic islet responses 
to glucose. Cell Metabolism, 24(3):389-401, 2016. 

7. Joan Camunas-Soler, Xiao-Qing Dai, Yan Hang, Austin Bautista, James Lyon, Kunimasa 
Suzuki, Seung K. Kim, Stephen R. Quake, and Patrick E. MacDonald. Patch-seq links single- 
cell transcriptomes to human islet dysfunction in diabetes. Cell Metabolism, 31(5):1017— 
1031.e4, May 2020. 

8. Eric A. Sobie. Parameter Sensitivity Analysis in Electrophysiological Models Using Multi- 
variable Regression. Biophysical Journal, 96(4):1264—1274, February 2009. 

9. Stefano Morotti and Eleonora Grandi. Logistic regression analysis of populations of electro- 
physiological models to assess proarrythmic risk. MethodsX, 4:25-34, January 2017. 

10. Michela Riz, Matthias Braun, and Morten Gram Pedersen. Mathematical modeling of het- 
erogeneous electrophysiological responses in human -cells. PLoS Computational Biology, 
10(1):e 1003389, January 2014. 

11. Matthias Braun, Reshma Ramracheya, Martin Bengtsson, Quan Zhang, Jovita Karanauskaite, 
Chris Partridge, Paul R Johnson, and Patrik Rorsman. Voltage-gated ion channels in hu- 
man pancreatic beta-cells: electrophysiological characterization and role in insulin secretion. 
Diabetes, 57(6): 1618-28, 2008. 

12. The MathWorks Inc. Global optimization toolbox version 4.8.1 (12023a), 2023. 

13. The MathWorks Inc. Parallel computing toolbox version 7.8 (12023a), 2023. 

14. Karoline Jeger and Aslak Tveito. Derivation of a Cell-Based Mathematical Model of Excitable 
Cells, pages 1-13. Springer International Publishing, 10 2020. 


5 Pancreatic Beta Cell Network 85 


15. 


16. 


17. 


18. 


19. 


20. 


21. 


22; 


23. 


24. 


25. 


26. 


27. 


28. 


29. 


30. 


31. 


32. 


Aslak Tveito, Karoline Jeger, Miroslav Kuchta, Kent-Andre Mardal, and Marie Rognes. A cell- 
based framework for numerical modeling of electrical conduction in cardiac tissue. Frontiers 
in Physics, 5, 10 2017. 

Karoline Jeger, Andrew Edwards, Wayne Giles, and Aslak Tveito. From millimeters to 
micrometers; re-introducing myocytes in models of cardiac electrophysiology. Frontiers in 
Physiology, 12, 10 2021. 

Ada Ellingsrud, Andreas Solbra, Gaute Einevoll, Geir Halnes, and Marie Rognes. Finite ele- 
ment simulation of ionic electrodiffusion in cellular geometries. Frontiers in Neuroinformatics, 
14:11, 03 2020. 

Ada Ellingsrud, Cécile Daversin-Catty, and Marie Rognes. A Cell-Based Model for Ionic 
Electrodiffusion in Excitable Tissue, pages 14-27. Springer International Publishing, 10 2020. 
Ashild Telle, James Trotter, Xing Cai, Henrik Finsberg, Miroslav Kuchta, Joakim Sundnes, 
and Samuel Wall. A cell-based framework for modeling cardiac mechanics. Biomechanics 
and Modeling in Mechanobiology, 22, 01 2023. 

Karoline Jæger and Aslak Tveito. Efficient, cell-based simulations of cardiac electrophysiology; 
the kirchhoff network model (knm). NPJ systems biology and applications, 9:25, 06 2023. 
Karoline Jeger and Aslak Tveito. The simplified kirchhoff network model (sknm); a cell-based 
reaction-diffusion model of excitable tissue. Unpublished, 2023. 

Gerardo J. Félix-Martínez, Aurelio N. Mata, and J. Rafael Godínez-Fernández. Reconstructing 
human pancreatic islet architectures using computational optimization. Islets, 12(6):1—13, 
2020. 

Gerardo J. Félix-Martínez and J. R. Godínez-Fernández. Comparative analysis of reconstructed 
architectures from mice and human islets. Islets, 14(1):23—35, 2021. 

Habo J. Jongsma and Ronald Wilders. Gap junctions in cardiovascular disease. Circulation 
Research, 86(12):1193-1197, 2000. 

J Takeda, M Gidh-Jain, L Z Xu, P Froguel, G Velho, M Vaxillaire, D Cohen, F Shimada, 
H Makino, and S Nishi. Structure/function studies of human beta-cell glucokinase. enzymatic 
properties of a sequence polymorphism, mutations associated with diabetes, and other site- 
directed mutants. J. Biol. Chem., 268(20):15200—-15204, July 1993. 

M Gidh-Jain, J Takeda, L Z Xu, A J Lange, N Vionnet, M Stoffel, P Froguel, G Velho, F Sun, 
and D Cohen. Glucokinase mutations associated with non-insulin-dependent (type 2) diabetes 
mellitus have decreased enzymatic activity: implications for structure/function relationships. 
Proceedings of the National Academy of Sciences, 90(5):1932-1936, March 1993. 

Peter M. Fernandes, James Kinkead, Iain McNae, Paul A.M. Michels, and Malcolm D. Walkin- 
shaw. Biochemical and transcript level differences between the three human phosphofructoki- 
nases show optimisation of each isoform for specific metabolic niches. Biochemical Journal, 
471(22):4425-4441, November 2020. 

Michael J. Mac Donald. Does glyceraldehyde enter pancreatic islet metabolism via both 
the triokinase and the glyceraldehyde phosphate dehydrogenase reactions? Archives of 
Biochemistry and Biophysics, 270(1):15—22, April 1989. 

Matthias Braun, Reshma Ramracheya, Martin Bengtsson, Quan Zhang, Jovita Karanauskaite, 
Chris Partridge, Paul R. Johnson, and Patrik Rorsman. Voltage-gated ion channels in hu- 
man pancreatic B-cells: Electrophysiological characterization and role in insulin secretion. 
Diabetes, 57(6):1618—1628, June 2008. 

R Bränström, Craig A Aspinwall, S Välimäki, C-G Óstensson, A Tibell, M Eckhard, H Brand- 
horst, BE Corkey, P-O Berggren, and O Larsson. Long-chain coa esters activate human pan- 
creatic beta-cell k atp channels: potential role in type 2 diabetes. Diabetologia, 47:277—283, 
2004. 

Barbara Rosati, Piero Marchetti, Olivia Crociani, Marzia Lecchi, Roberto Lupi, Annarosa Ar- 
cangeli, Massimo Olivotto, and Enzo Wanke. Glucose-and arginine-induced insulin secretion 
by human pancreatic B-cells: the role of HERG K+ channels in firing and release. The FASEB 
Journal, 14(15):2601—2610, 2000. 

Stefan Meier, Adaia Grundland, Dobromir Dobrev, Paul G. A. Volders, and Jordi Heijman. In 
silico analysis of the dynamic regulation of cardiac electrophysiology by kv11.1 ion-channel 
trafficking. The Journal of Physiology, 601(13):2711—2731, February 2023. 


86 Pancreatic Beta Cell Network 


33. Brodie A. J. Lawson, Christopher C. Drovandi, Nicole Cusimano, Pamela Burrage, Blanca 
Rodriguez, and Kevin Burrage. Unlocking data sets by calibrating populations of models to 
data density: A study in atrial electrophysiology. Science Advances, 4(1), January 2018. 

34. Trine Krogh-Madsen, Eric A. Sobie, and David J. Christini. Improving cardiomyocyte model 
fidelity and utility via dynamic electrophysiology protocols and optimization algorithms. The 
Journal of Physiology, 594(9):2525—2536, 2016. 

35. Trine Krogh-Madsen, Anna F. Jacobson, Francis A. Ortega, and David J. Christini. Global 
Optimization of Ventricular Myocyte Model to Multi- Variable Objective Improves Predictions 
of Drug-Induced Torsades de Pointes. Frontiers in Physiology, 8, 2017. 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, 
distribution and reproduction in any medium or format, as long as you give appropriate credit to the 
original author(s) and the source, provide a link to the Creative Commons license and indicate if 
changes were made. 

The images or other third party material in this chapter are included in the chapter's Creative 
Commons license, unless indicated otherwise in a credit line to the material. If material is not included 
in the chapter's Creative Commons license and your intended use is not permitted by statutory 
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright 
holder. 


eo 
Chapter 6 E 
The Impact of Mechano-Electric Feedback on 
Drug- and Stretch- Induced Arrhythmia Using a 
Computational Model of Cardiac 
Electromechanics 


Anthony Asencio, Melania Buonocunto, Matthew W Ellis, Karl Munthe, Kyle T 
Stark, Joakim Sundnes, Henrik Finsberg, and Hermenegild J Arevalo 


Abstract Mechano-electric feedback (MEF) is thought to be an important factor 
in the increased arrhythmic risk observed clinically in heart failure and chronic 
infarct patients. Here, we utilize pathologic stretch as a parameter for investigating 
the sensitivity of cardiac electrophysiology and mechanics to MEF. To simulate 
the effects of stretch, we incorporated stretch-activated ion channels (SACs) into 
a coupled electro-mechanics model based on the well-established O'Hara-Rudy 
(ORd) (electrophysiology) Land (mechanics) models. We investigated the effect 
different degrees of stretch had upon electrophysiological parameters such as the 
action potential duration and calcium transient, as well as functional parameters 
like force production. We further determined the sensitivity of cardiac cells under 
stretch to various pro- and anti-arrhythmic drugs to better inform on drug risk 
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classifications and indications. Our results indicate that pathologic stretch may exert 
a severe pro-arrhythmic effect on cardiac cells, which markedly exacerbates the risk 
profile of pro-arrhythmic drugs, but one which was mitigated through the action of 
anti-arrhythmic compounds. 


6.1 Introduction 


Patients with heart failure have increased risk of arrhythmia, which may be ex- 
acerbated by some medications used in treatment. The Comprehensive in vitro 
Pro-Arrhythmia Assay (CiPA) involves computational and experimental modeling 
alongside clinical observations to form a risk stratification profile for clinically 
prescribed medications for heart failure patients. In the present study we sought 
to further inform these risk profiles by incorporating mechano-electric feedback 
(MEF) into the in silico modeling of drug effects, using a coupled model based on 
the O'Hara-Rudy model (ORd) for human ventricular electrophysiology [1] and the 
model by Land et al. for myocyte mechanics [2]. The coupling of the two models 
occurs primarily through intracellular calcium, which regulates tension development 
by binding to Troponin C, see, for instance, [2] for details. This binding is tension 
dependent, and thereby creates an MEF effect where the mechanical state of the 
myocyte affects the concentration of free calcium and thereby the electrical activity 
of the cell. 

MEF is thought to be a significant contributor to spontaneous arrhythmia de- 
velopment in heart failure patients [3], but has yet to be deeply investigated in 
computational studies coupling cardiac electrophysiology and mechanics. In addi- 
tion to the tension-dependent calcium buffering, known mechanisms of MEF include 
stretch-activated ion channels (SACs) [4]. These are channels in the cell membrane 
that open in response to mechanical stimuli (i.e., cell stretch), and thereby allow 
ions to flow in and out of the cell. Three SACs are typically recognized in car- 
diac tissue; non-selective cation channels (SACns) primarily conducting sodium and 
potassium (and, to a lesser degree, calcium), potassium-selective channels (SACK), 
and calcium-selective channels (SACCa). We incorporated established and physi- 
ologically representative mathematical formulations for these three SACs into the 
coupled electro-mechanics model. 

Furthermore, we chose various CiPA compounds across the risk profile, includ- 
ing Dofetilide and Bepridil (high risk), Domperidone and Cisapride (intermediate 
risk), and Ranolazine and Diltiazem (low risk). These drugs were chosen due to their 
role in blocking either rectifying potassium currents, which contributes to elongated 
action potentials (AP) and possible early after depolarization events, calcium cur- 
rents, a critical linkage mechanism in MEF, or both. We sought to assess the effect 
of pathological stretch via SAC activation on electrophysiology and mechanics at 
baseline, and in combination with drug activity to investigate any altered sensitivity 
to drug action. 
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6.2 Models and Methods 


In this section we first describe the single cell electro-mechanics model, including 
MEF and the modifications to incorporate drugs, as well as the experimental protocol. 
The coupled electro-mechanics model was based on the published version of the ORd 
and Land et al models [1, 2], and we here focus on the coupling of the models in terms 
of calcium binding, force development, and SACs. For a complete specification of 
the models we refer to the original publications. 


6.2.1 Stretch-Activated Channels 


We incorporated three SACs currents in the coupled ORd-Land model: Js 4cns, 
IsAcKo and JsAccap. The implementation of the first two currents was based on 
the formulation by Niederer and Smith (2007) [5], while their calibration and the 
stretch dependence of the Js 4cns expression were based on Buonocunto et al. (2023) 
[6]. Finally, we formulated the last current deriving it from the expression of the 
background calcium current of the ORd model. The mechanical input to the SACs is 
the extension ratio A, which is the ratio of the current sarcomere length (SL) to the 
slack length. The SACs are formulated so that they are closed and non-conducting 
for A < 1 (contraction), and open for A > 1 (stretch). 

The stretch-activated /sAc4; current was defined as the sum of a sodium and a 
potassium component: 


IsACns = IsACnsyq + ISACnsy: (6.1) 

with: 
Ts aCnsna =F : 8ns* YSL,ans * (Vin a Ena); (6.2) 
IsACnsk = gns " YSL,K ^ (Vin — Ex). (6.3) 


Here, gns is the maximal conductance of the SAC,,, ENa and Ex are the reversal 
potentials of sodium and potassium, and YsŁ,ns the sigmoidal stretch dependence of 
the conductance, defined as 


fas (A — 1 
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Furthermore, the stretch-activated Igo current was formulated as 


SL,Ko ` (V — Ex) 
IsACKo = 8Ko* 7 2 Tov , (6.4) 
l+e 5 


with gxo being the maximal conductance of the SACK o, and YsL,Ko the stretch 
dependence of the conductance defined by YsL,Ko = £o (A — 1) + 0.7. Finally, we 
formulated the the calcium mechanosensitive current as 
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Isaccap = 4'Gcap:x:(-0.341-[Ca?*], - [Ca?*];-e RT) — 77 — ——, (6.5) 


with GCap being the channel conductance, F, R, and T conventional thermody- 
namic constants, [Ca?*]; and [Ca?*], the intracellular and extracellular calcium 
concentrations. In addition, both Is 4cys and /sAcko were divided by the membrane 
capacitance C,, which was expressed as a Hill-curve that increases with increasing 
stretch: 

(A - p* 
0.0354 + (a - 1)* 


as experimentally observed by [7] and parameterised by [6]. 


Cm = 1.534 x 107+ , (6.6) 


6.2.2 Mechanical Modeling 


The mechanical aspect of our model was addressed by the incorporation of the 
Land model of cardiac contraction [2]. Here, the total force is expressed as a sum 
of different passive and active contributions [8]. The active force results from the 
crossbridge cycle, in which ATP-driven myosin translates along filamentous actin. 
The resulting force is expressed as the product of the mean cross-bridge distortion 
and cross-bridge stiffness, assuming that these behave as elastic springs. Passive 
forces include an elastic component coming mainly from titin and collagen, and a 
viscous component. Force balance requires that the net force on a myocyte is zero: 


Fa + Fp + Fse + Fy + Fore =0, (6.7) 


where F4 is the active force, F is the passive elastic force, Fse = Kse(A — Aset) is 
a series elastic force, for which the spring stiffness Kse and set stretch Aset can be 
varied to mimic different experimental settings. Furthermore, F, is the viscous force 
and Fpre is a constant preload applied to the myocyte. 

The formulation of the active force is taken from [2], and given by 


Tye r 
Fa = WAZE + &)8 + GW). (6.8) 


Here, Tep is the tetanic (maximum) force at the resting length, A = 1 [8, 2], rs is 
the steady state duty ratio, while A(4) is a phenomenological function that accounts 
for the force-length relationship. Finally, the expression ((1-- ¢s)S - Z4,W) describes 
how the dynamic force varies with the crossbridge cycling, with Z, and Zy,, being 
the distortion of the crossbridges in the strongly (S) and weakly (W) bound states, 
respectively. 

The sensitivity of active force to stretch is increased by the cooperativity of 
tropomyosin along actin. The binding of calcium to a single troponin site increases 
the probability of calcium binding to the neighboring sites, thereby increasing the 
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myosin binding probability and, as a result, the active force. For a more detailed 
description and complete formulation of the equations, we refer to Land, et al. [2]. 


6.2.3 Modeling the Effect of Drugs 


To model the activity of various drugs on ion channel conductance, we used the 
following equation: 


SF 


_ Bdrug _ | i (parece a 


g ICso 


Here g is a measure of ion channel conductance, EFTPC is the effective free 
therapeutic plasma concentration of each drug, ICs is the inhibitory concentration 
of the drug necessary to bring each respective channel activity down to 50%, and A is 
the Hill coefficient characterizing cooperativity of the drug and channel interaction. 
EFTPC, ICso, and h for each drug and channel used in this study were based on 
previously published experimental and clinical results [9, 10, 11]. SF indicates the 
scaling factor for each respective channel, which could be used as a parameter in the 
coupled model to simulate drug action. 

As an example, Table 6.1 summarizes the scaling factors for the various ion 
channels for the compound Ranolazine. 


Table 6.1: Ion channel scaling factors (SF) for various plasma concentrations of 
Ranolazine 


Channel SF 1x EFTPC 2x EFTPC 4x EFTPC 
Ikr 0.810 0.681 0.516 
Iya 0.663 0.512 0.360 
INyaL 0.800 0.668 0.503 
Icar 0.975 0.954 0.918 
Iks 0.994 0.989 0.978 
Iki 1.000 1.000 1.000 
lho 1.000 1.000 1.000 


6.2.4 Experimental Protocol 


After pre-pacing our simulations for 50 beats to reach steady state for all the 
electrophysiological signals, we ran isometric test simulations using the govern- 
ing equations presented above and Kse = 10°. The stretch was kept constant (at 
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A= 0.9,0.95, 1, 1.05, 1.1, 1.2) for the duration of the cycle length (1000 ms) while 
the electromechanics signals were evaluated. Furthermore, to observe the effect of 
quick stretch release on the electromechanics signals, we performed a stretch/release 
test. It consisted of holding stretch constant at a set value (A = 1.00, 1.05, 1.1) for 600 
ms, and then suddenly releasing 2% of the applied stretch and keeping it constant 
for the remaining cycle length (400 ms). 


6.3 Results 
6.3.1 Isometric Behavior 


Application of isometric stretch and contraction in presence of SACs did not provoke 
any changes when 2 < 1, while it caused prolongation of the AP with A > 1 
(figure 6.1 A). Additionally, 4 = 1.20 caused a complete loss of AP morphology and 
prevented repolarization. Similarly, the intracellular Ca?* peak value increased with 
all A values, although it provoked larger increase with A > 1 (figure 6.1B). However, 
A = 12 completely flattened the signal. Similar behaviour was observed for the Ca?* 
bound to Troponin (data not shown). Furthermore, the active force signal was found 
to increase positively with both expansion (stretch > 1) and compression (stretch < 
1) in figure 6.2A in absence of SACs. The addition of current from SACs (figure 
6.2B) prolonged the duration of the active force signal and accentuated the dome 
shape for values of stretch > 1. Stretch values less than or equal to one were not 
effected by the presence of SACs. However stretch > 1 increased active force relative 
to the no-SACs case. 
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Fig. 6.1: Membrane Potential A) and Intracellular [Ca?*] B) over a range of isometric 
stretches in presence of SACs. 
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Fig. 6.2: Active forces over a range of isometric stretches A) without SACs and B) 
with SACs. 


6.3.2 Stretch/Release Test 


We performed the stretch/release test and evaluated the changes to the electrome- 
chanics signals at the moment of release and in presence of SACs. The membrane 
voltage underwent a small (< 5%) sudden decrease in amplitude (figure 6.3A), while 
the active force showed a more pronounced drop (maximum 25%) and a subsequent 
recovery of the signal with decreased amplitude (figure 6.3B) with 4 > 1. 
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Fig. 6.3: Effects of stretch/release test on A) Membrane potential and B) Active force 
in presence of SACs. 
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6.3.3 Drug Effects in Isometric Conditions 


Using (6.9) with reference values for drug concentration and activity [9, 10, 11], 
we assessed the effect of each drug on the AP and force output. The results are 
summarized in Table 6.2, and a representative plot is shown in figure 6.4. 
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Fig. 6.4: Effect of various CiPA drug compounds on the cardiac action potential. 


Table 6.2: Effect of CiPA drugs on the action potential duration at 90% repolarization 
AP Doo and maximum force 


Drug APDop (ms) Maximum Force 
(kPa) 
Control 263.7 20.6 
Dofetilide 317.9 22.8 
Bepridil 339.2 21.5 
Domperidone 296.2 22.5 
Cisapride 312.8 22.6 
Ranolazine 286.1 28.0 


Diltiazem 236.9 10.3 
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Of the drugs assessed, Bepridil displayed the most pronounced effects on the AP 
relative to the control state, causing the largest increase of APDog at the EFTPC. To 
further classify the pro-arrhythmic effect of Bepridil, we ran additional simulations 
at twice and four times the EFTPC of the drug, and observed worsening on the 
above outcomes, including longer APDoo and lower maximum force values with 
higher doses, as well as the appearance of early after depolarizations indicative of 
arrhythmia risk. 


6.3.4 Combining Drugs with Stretch 


In order to evaluate the effect of MEF more comprehensively, we next combined our 
modeling of pathologic stretch together with drug implementations. Based on known 
drug activity and our results at the EFTPC of each simulated drug, we focused our 
efforts on Bepridil (a high-risk CiPA pro-arrhythmic compound which elongated 
APDo, versus control) and Diltiazem (a low-risk CiPA anti-arrhythmic compound 
which shortened the APDog and greatly reduced maximum force versus control) to 
characterize the extent to which stretch sensitizes cardiomyocytes to drug action. As 
shown in Figure 6.5, we found that the introduction of 5% pathologic stretch served 
to elongate the APDoo at 1x EFTPC Bepridil by over 200 ms versus the unstretched 
condition, a duration twice that of the control condition, while exacerbating the 
appearance of early after depolarization events. On the other hand, the presence of 
Diltiazem greatly mitigated the deleterious effects of pathologic stretch, with the 1x 
EFTPC and 5% pathologic stretch maintaining an action potential of similar duration 
and appearance to the control condition. 
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Fig. 6.5: Effect of 2X plasma concentration of the pro-arrhythmic Bepridil and anti- 
arrhythmic Diltiazem on the cardiac action potential when combined with pathologic 
stretch. 
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6.4 Conclusion 


We have developed a coupled model of mechano-electric feedback in human car- 
diomyocytes, by incorporating mechano-electric feedback in the form of SACs in the 
ORd-Landcoupled electromechanics model. The SACs were formulated so that they 
would have no effect on sarcomeric force production or AP morphology when stretch 
was less than or equal to one, which is consistent with previous studies. However, 
under conditions of pathological stretch, the SACs were shown to amplify the force- 
stretch relation and the increased Ca?* binding with stretch, which is expected since 
the SACs give rise to an increase in Ca?* influx. Following this finding, we incor- 
porated non-physiologic stretch into our drug implementation scheme. Focusing on 
the most pro-arrhythmic compound assessed, Bepridil, and the most anti-arrhythmic 
drug assessed, Diltiazem, we were able to characterize the interaction between drug 
activity and pathologic stretch. Notably, we found that even 5% non-physiologic 
stretch was able to markedly exacerbate the pro-arrhythmic status of Bepridil, elon- 
gating the APDoọ significantly, producing AP signatures more similar to many times 
the effective free therapeutic plasma concentration of the drug. This was in contrastto 
the anti-arrhythmic calcium channel blocker Diltiazem, which served to mitigate the 
pathological activity of non-physiological stretching, maintaining an AP signature 
highly similar to the unstretched control condition. The ORd-Land model with SACs 
was also straight forwardly implemented in a 3D finite element method (FEM) solver 
based on FEniCS! and used to compute the contraction and relaxation of a small 
rectangular slab of heart muscle tissue. Preliminary results indicate that the tissue 
reacts qualitatively similar to the cell model, but further investigation is needed to 
quantify how the 3D FEM model changes when SACS are included and compare the 
simulation results to experimental results. As mechano-electric feedback remains 
an understudied consideration in arrhythmia risk stratification models, we believe 
our findings here provide a more comprehensive model for cardiac disease modeling 
and better inform drug risk classifications and indications, thereby improving clinical 
care for patients with heart disease. 


References 


1. Thomas O'Hara, László Virág, Andras Varró, and Yoram Rudy. Simulation of the undiseased 
human cardiac ventricular action potential: model formulation and experimental validation. 
PLoS computational biology, 7(5):e1002061, 2011. 

2. Sander Land, So-Jin Park-Holohan, Nicolas P Smith, Cristobal G Dos Remedios, Jonathan C 
Kentish, and Steven A Niederer. A model of cardiac contraction based on novel measurements 
of tension development in human cardiomyocytes. J. Mol. Cell. Cardiol., 106:68—83, May 
2017. 

3. T Alexander Quinn and Peter Kohl. Cardiac mechano-electric coupling: acute effects of 
mechanical stimulation on heart rate and rhythm. Physiological reviews, 101(1):37—92, 2021. 


! https://fenicsproject.org/ 


6 Cardiac Mechano-Electric Feedback 97 


4. Rémi Peyronnet, Jeanne M Nerbonne, and Peter Kohl. Cardiac mechano-gated ion channels 
and arrhythmias. Circulation research, 118(2):311—329, 2016. 

5. Steven A Niederer and Nicolas P Smith. A mathematical model of the slow force response to 
stretch in rat ventricular myocytes. Biophysical journal, 92(11):4030-4044, 2007. 

6. Melania Buonocunto, Aurore Lyon, Tammo Delhaas, Jordi Heijman, and Joost Lumens. Elec- 
trophysiological effects of stretch-activated ion channels: a systematic computational charac- 
terization. The Journal of Physiology, 2023. 

7. Bernardo L de Oliveira, Emily R Pfeiffer, Joakim Sundnes, Samuel T Wall, and Andrew D 
McCulloch. Increased cell membrane capacitance is the dominant mechanism of stretch- 
dependent conduction slowing in the rabbit heart: a computational study. Cellular and molecular 
bioengineering, 8:237—246, 2015. 

8. A. Raaum. Computational study of the impact of parameter variability and drugs on healthy 
and failing cardiac cells: A population of models approach to characterize drug effects on cell 
electrophysiology and mechanics. University of Oslo, Spring 2023. 

9. William J. Crumb, Jose Vicente, Lars Johannesen, and David G. Strauss. An evaluation of 
30 clinical drugs against the comprehensive in vitro proarrhythmia assay (cipa) proposed ion 
channel panel. Journal of Pharmacological and Toxicological Methods, 81:251—262, 2016. 
Focused Issue on Safety Pharmacology. 

10. Zhihua Li, Sara Dutta, Jiansong Sheng, Phu N. Tran, Wendy Wu, Kelly Chang, Thembi Mdluli, 
David G. Strauss, and Thomas Colatsky. Improving the in silico assessment of proarrhythmia 
risk by combining herg (human ether-à-go-go-related gene) channel-drug binding kinetics and 
multichannel pharmacology. Circulation: Arrhythmia and Electrophysiology, 10(2):e004628, 
2017. 

11. Jordi Llopis-Lorente, Beatriz Trenor, and Javier Saiz. Considering population variability of 
electrophysiological models improves the in silico assessment of drug-induced torsadogenic 
risk. Computer Methods and Programs in Biomedicine, 221:106934, 2022. 


Open Access This chapter is licensed under the terms of the Creative Commons Attribution 4.0 
International License (http://creativecommons.org/licenses/by/4.0/), which permits use, sharing, adaptation, 
distribution and reproduction in any medium or format, as long as you give appropriate credit to the 
original author(s) and the source, provide a link to the Creative Commons license and indicate if 
changes were made. 

The images or other third party material in this chapter are included in the chapter's Creative 
Commons license, unless indicated otherwise in a credit line to the material. If material is not included 
in the chapter's Creative Commons license and your intended use is not permitted by statutory 
regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright 
holder. 


S 
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Impact of Modeling Assumptions on 
Hemodynamic Stresses in Predicting Cerebral 
Aneurysm Rupture Status 


Guillermo L Nozaleda, Sofia Poloni, Luca Soliveri, Kristian Valen-Sendstad 


Abstract Approximately 3% of the population is estimated to have cerebral 
aneurysms, which are the leading cause of subarachnoid haemorrhage. Convinc- 
ing evidences suggest that wall shear stresses (WSS) play a role in vessel remodeling 
and in the development of vascular diseases. Since WSS cannot be directly measured, 
researchers have resorted to using medical images available in routine clinical prac- 
tice to simulate computational fluid dynamics (CFD) and investigate patient-specific 
vascular conditions. They retrospectively analyse the correlation between WSS and 
disease outcomes to find potential clinical tools for future use. However, most of 
these models are based on assumptions that introduce variability and error. In this 
work we investigated the effects of a non-Newtonian viscosity model and inflow 
uncertainty on the prediction of commonly computed hemodynamic metrics. Our 
results show a substantial influence of the non-Newtonian model and blood flow rate 
on CFD outcomes, highlighting the need of incorporating non-Newtonian rheology 
and patient-specific blood flow measurements in CFD simulations. 
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7.1 Introduction 


Roughly 3% of human population is estimated to be affected by cerebral aneurysms, 
which stand as the primary trigger for subarachnoid hemorrhages [1]. Aneurysm 
formation involves a complex interplay between biological processes within the blood 
vessel wall and hemodynamic factors [2]. The underlying hypothesis suggests that 
abnormal WSS triggers an inflammatory reaction, impacting arterial wall remodeling 
and significantly contributing to aneurysm initiation and progression [3]. If the 
arterial wall fails to withstand these forces, the aneurysm may rupture, resulting in 
high risk of mortality [4]. 

In the management of an aneurysm, evaluating risk of rupture is thus of utmost 
importance for clinicians [5]. Hemodynamic factors, especially WSS measurements, 
serve as common rupture indicators. However, direct WSS measurement from clini- 
cal observations is not feasible [6]. As an alternative, researchers employ real geome- 
tries from routine clinical practice to perform computational fluid dynamics (CFD) 
simulations and study patient-specific vascular conditions [7, 8, 9, 10, 11, 12]. Nev- 
ertheless, the reliability of CFD in the context of aneurysm management remains 
a subject of ongoing debate. In this context, a recent review by Steinman [6] ex- 
plores each stage of the patient-specific CFD pipeline, emphasizing the influence of 
sources of variability and error. These sources include imaging techniques, blood 
rheology and artery wall characterization, prescribed boundary conditions, segmen- 
tation methods, and variations in solver numerics. While some of this variability can 
be mitigated through appropriate modeling, some assumptions are more critical to 
overcome, primarily due to the lack of patient-specific data. 

Blood behaves as a non-Newtonian fluid and artery walls undergo cyclic deforma- 
tions. Nevertheless, Steinman [6] argues that the potential benefits of including these 
characteristics in CFD simulations might be outweighed by the computational ef- 
fort required. Consequently, many studies prefer simplified assumptions of constant 
viscosity and rigid walls, reducing computational cost and modeling complexity. 
Despite some efforts to investigate the effects of these simplifications [13, 14, 15], 
the debate surrounding proper modeling techniques remains open. When it comes to 
boundary conditions, it is important to establish reasonable outflow rates by employ- 
ing proper scaling laws [16]. As for inflow rates, the generalized absence of patient- 
specific data makes it necessary to rely on estimations. This becomes a primary 
limitation, particularly since WSS magnitude and distribution are highly responsive 
to variations in these prescribed values [13, 17]. Additionally, careful consideration 
should be given to the locations of model truncation and diameter measurements to 
ensure accurate results [6]. Regarding solver numerics, Valen-Sendstad [18] demon- 
strated how commonly used schemes, meshes, and time steps in aneurysm CFD may 
dampen flow instabilities and alter stress distribution compared to higher-order ap- 
proaches adopting fine spatial and temporal discretisations. Khan et al. [19], on the 
other hand, warned against relying on default schemes from commercial software. 
Furthermore, many studies assumed laminar flow and employed dissipative schemes 
for stability, which can potentially affect accuracy. Therefore, achieving a reliable 
numerical solution entails advanced schemes, avoiding laminar assumptions, ensur- 
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ing minimal dissipation, and fine discretisations. Our work considers these aspects 
for improved precision and robustness. 

In the literature, most studies tend to focus on assessing the impact of a single 
factor on WSS. This often involves comparing works carried out with different 
CFD settings and methodologies, which may lead to ambiguous and misleading 
conclusions due to their reliance on particular numerical aspects. To address this 
limitation, our present study takes a more comprehensive approach by conducting 
an analysis that considers various modeling assumptions. The primary objective is 
to evaluate how these different conditions can influence the CFD outcomes, thereby 
providing further insights into the sources of variability and error. Specifically, the 
current work compares the impact of the controversial non-Newtonian model with 
that of inflow rates, which is universally considered to be one of the most influential 
between the modeling assumptions. 


7.2 Methods 
7.2.1 Study Population and Design 


Five aneurysms geometries were selected form the open source AneuriskWeb repos- 
itory [20]. These cases all involve internal carotid artery (ICA) aneurysms to ensure 
analysis consistency. For each geometry, we conducted four simulations: i) standard 
simulation (VO) with an standard inflow and constant viscosity, ii) non-Newtonian 
viscosity model, iii) 25% reduction in blood flow, and iv) 25% increase in blood 
flow. The conducted simulations, based on the standard deviations reported by Hoi 
et al. |21], were solved using the high performance Saga HPC cluster, a compute 
node with 20 Intel Xeon-Gold 6138 2.0 GHz cores. The typical computation time 
per simulation averaged 12 hours. 


7.2.2 Computational Fluid Dynamics 


Consider an incompressible fluid of kinematic viscosity v and density p, in the 
absence of body forces the problem reduces to the integration of 


P eu V)u = vV?u - Vp, (7.1) 


V.u-0, (1.2) 


where u = (u1, u2, ua) and p are the velocity and reduced pressure of the fluid. 

The Navier-Stokes equations (7.1)-(7.2) were numerically solved using Oasis, 
an open-source and validated CFD solver [22]. Oasis is based on the Finite Ele- 
ment Method implemented on the FEniCS computing platform [23]. It utilizes a 
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segregated, space/time centered, incremental pressure correction scheme to achieve 
accurate solutions while minimizing numerical dispersion and diffusion [24]. Below, 
we provide a description of the spatial and temporal discretisations used in Oasis. 
For more detailed information about the solver and its numerical implementation, 
readers are referred to the original manuscript by Mortensen et al. [22]. 


7.2.2.1 Fractional Step Algorithm 


In Oasis, a space/time catered fractional step method is used to advance in time. Let 
[0,7] represent a cardiac cycle, which is divided into uniform intervals of length 
At = t — t! with n € Z. Specifically, we set 7 = 0.951 s and At = 0.0951 
ms, resulting in 10,000 timesteps per cycle. Simulations spanned 3 cardiac cycles. 
Denoting u% as the k velocity component at time /", the algorithm looks as follows 


ul _ yr! 
Step 1: | EE um +B? WR -Vyp' for k-(L2,), (73) 
1 
Step 2: Vy = d ul, (7.4) 
ug = u; 
Step | APO Ve fo k={1,2,3}, (7.5) 


where y = p”~!/? — p* denotes a pressure correction and p* is a tentative pressure. 
In the first step (7.3), the viscous terms are handled using a Crank-Nicolson scheme, 
such that ij = 0.5 (ui, + ur. and the convective terms are incorporated in BY 2 


For the latter computation, Oasis provides two options, 


nO? it- Vin. (1.6) 


p". 3,1 Vu! — le» Vu? and B 
The first approach in (7.6) consists of a fully explicit Adams-Bashforth discretization. 
In contrast, the second scheme is implicit, with 4 = 1.5u”7! — 0.5u"-? and ak = 
0.5 (uj, + ut), as given above. 

Our objective is to calculate the updated values of u” and p"-!/? for each time 
step. To accomplish this, we start by solving (7.3) to determine the tentative velocity 
ul, which is non-solenoidal (ie., V - u” + 0). The velocity u! is then used to 
determine p"-!/? upon resolution of Poisson’s equation (7.4). Finally, we calculate 


from (7.5) the updated divergence-free velocity field u". 
7.2.2.2 Variational Formulation 
For the spatial discretization, the finite element method is used to discretise (7.3)— 


(7.5) on the bounded domain Q c R2, with boundary ðQ. All simulations em- 
ployed second-degree Lagrange polynomials for velocity and first-degree for pres- 
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sure. This scheme guarantees second-order accuracy in the L2 norm, leading to an 
error reduction rate of at least O(h?) as the mesh size h is refined. Meshes with 
210, 000 — 240, 000 elements were obtained using VaMPy, which is equivalent to 
1.68 — 1.92 million linear elements [25]. We define the trial and test spaces 


V-(veH'(Q):v2uoonóQ], V= {ve H'(Q):v=00ndQ}, (7.7) 


where H! (Q) is the Sobolev space, which contains functions v satisfying that v? 
and |Vv|? have finite integrals over Q. Consider v and q to be the test functions for 
the velocity and pressure, respectively. The variational formulation of (7.3)-(7.5) 
can be obtained multiplying each equation by the corresponding test function and 
integrating along the domain, thereby yielding 


ut » ul n-1/2 a = x * 

—— + B, v+vVi,: Vvdx = yV,üyvds — | Vep vdx, (7.8) 
Q At ðQ Q 
V-ul 
Vy: Vqdx - Ving ds = q dx, (7.9) 
Q aQ o At 
ug -uk 
f —— — 7 y dx =- vernis, (7.10) 
o At Q 


after use is made of integration by parts. Above, V,, represents the gradient in the 
outward normal direction to the boundary 0Q. 


7.2.2.3 Boundary Conditions and Viscosity Models 


In the absence of patient-specific measurements, an average blood flow of 4.05 cm? /s 
was applied at the inlet, referencing [21]. The inlet boundary conditions were estab- 
lished using the generalized blood flow curve from Oasis. Outlet pressure conditions 
were adjusted to achieve a blood flow split based on outlet cross-sectional areas. 
Blood's viscosity changes in response to shear stress, a phenomenon known 
as non-Newtonian behavior. Despite this knowledge, many studies on intracranial 
aneurysms simplify blood as Newtonian fluid [10, 11, 12], assuming constant viscos- 
ity. This work aims to contrast the Newtonian approximation with a more accurate 
representation of the fluid rheology. In particular, the modified Cross model [26, 27] 


LA) = Hoo _ 1 (711) 


Ho- Hoo — [1+ (ay)"]* 


is adopted here to characterize blood’s non-Newtonian behavior. In the above expres- 
sion y is the shear rate, {fo = 3.372 mPa - s and uo = 3.372 mPa - s are the infinite- 
shear and zero-shear viscosities, respectively, whereas 4 = 3.736 s, m = 2.406, and 
a = 0.34 represent fitting parameters. The kinematic viscosity, used in the above 
formulation, can be obtained from (7.11) as v = u/p. 
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7.2.3 Post-Processing 


The post-processing of the results was performed exploiting VaMPy [25] func- 
tionalities and Paraview software. The post-processing was performed on the third 
cycle, to avoid artifacts due to initialization, and focusing on the aneurysm sac only. 
Commonly used rupture risk metrics, specifically time-averaged wall shear stresses 
(TAWSS) and oscillatory shear index (OSD, were calculated as follows 


T 
1 |f, vwdtl 


1 T 
TAWSS = = I twdt and OSI- --H (7.12) 
T Jo 2 fy lrwlat 


Furthermore, VaMPy was adjusted to handle non-Newtonian effects during post- 
processing, including the storage of viscosity values from simulations for stress 
calculations. 


7.3 Results 


Figure 7.1 illustrates the qualitative velocity magnitude isosurfaces at the systolic 
peak for the four cases under analysis: the standard simulation (VO), the non- 
Newtonian model, and inflow rates altered by +25%. In the standard simulation, 
distinct low-velocity recirculation zones and a complex flow phenotype are evident 
within the aneurysm sac. When incorporating non-Newtonian effects, the velocity 
gradient in the same region decreases due to higher viscosity values associated with 
low shear rates. In scenarios where the inflow rate is decreased or increased by 25%, 
there is a corresponding reduction or augmentation of velocity within the aneurysmal 
area. 

Examining the distribution of TAWSS in Figure 7.2A, it is evident that TAWSS 
registers higher values at the neck of the aneurysm and exhibits lower values within 
the sac, where velocity recirculation is prominent. This pattern remains consistent 
for all models considered here. Furthermore, Figure 7.2B quantifies the effect of 
modeling assumption on the TAWSS within the aneurysm sac with respect to the VO 
simulation. In the case of the non-Newtonian assumption, stresses decrease due to 
heightened viscosity. This effect mirrors the consequences of reduced inflow, where 
it is intuitive that lower flow results in decreased TAWSS. Likewise, an increase in 
inflow leads to a rise in TAWSS within the aneurysm sac, as expected. To conclude 
this part of the analysis, Table 7.1 outlines the influence of different modeling 
assumptions on TAWSS, expressed as a percentage alteration with respect to the VO 
simulation. The numbers here reported solidify the observation that all cases exert 
a substantial impact on TAWSS. Notably, TAWSS displays heightened sensitivity to 
variations in flow beyond what would be expected from a linear relationship. The 
influence of the non-Newtonian assumption on stress levels is particularly significant, 
resulting comparable to a 25% decrease in inflow. 
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Fig. 7.2: a) Maps of TAWSS distribution. b) Effect of modeling assumptions on 
TAWSS with respect to the VO simulation. Q: blood flow. TAWSS: time average wall 
shear stress. 
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Table 7.1: Impact of modeling assumptions: percentage variation of TAWSS with 
respect to VO simulation. Q: blood flow. TAWSS: time average wall shear stress. 


TAWSS 
Non-Newtonian Q-2596 Q 42596 
C0001 -37,62 -38,13 44,01 
C0002 -41,76 -38,69 43,17 
C0005 -32,50 -33,10 32,98 
C0006 -37,29 -37,40 40,01 
C0067 -39,58 -35,20 39,15 


Attention is given in Figure 7.3A to the OSI distribution for all geometries and 
models, highlighting regions within the aneurysmal sac with elevated OSI values, 
indicative of disturbed flow. From the quantitative results (Figure 7.3B and Table 
7.2), it follows that OSI variations are more complex and do not correlate with the 
variations of inflow. Instead, the pronounced variability suggests that OSI is more 
sensitive to model geometries. Concerning the impact of non-Newtonian effects, the 
results underscore a substantial influence on OSI magnitude. 


a) vo Non-Newtonian 


C0002 


C0006 


C0067 


BF -25% 


BF+25% b) 


ovo 
© Non-newtonian 
= 0-255 W coos 
+ Q425% 


ago 
FF To001 
= 


oos <> 


o 
002 + = 
O co067 
= 
001 o 
o 
0 
0,01 0,02 0,03 0,04 0,05 0,06. 0,07 0,08 
OSI 
0 osi 0.5 
— -— 


Fig. 7.3: a) Maps of OSI distribution. b) Effect of modeling assumptions on OSI 
with respect to the VO simulation. Q: blood flow. OSI: oscillatory shear index. 
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Table 7.2: Impact of modeling assumptions: percentage variation of OSI with respect 
to VO simulation. Q: blood flow. OSI: oscillatory shear index. 


OSI 
Non-Newtonian Q- 2596 Q+25% 
C0001 -41,29 -8,53 4,64 
C0002 -25,40 -6,15 0,11 
C0005 -64,08 -32,46 23,23 
C0006 -23,09 1,23 -0,08 
C0067 -56,03 -20,85 22,42 


7.4 Conclusion 


7.4.1 Impact 


Medical images available from clinical practice are extensively utilized for simu- 
lating CFD and exploring patient-specific vascular conditions, with the final aim 
of predicting rupture in cerebral aneurysms. However, most computational models 
are patient-specific in terms of their lumen geometry only and rely on numerous 
assumptions, which introduce a wide range of uncertainty and error. 

In the current study we have performed CFD simulations under different mod- 
eling assumptions on 5 ICA aneurysms geometries. Specifically, our work aims at 
comparing the impact of the controversial non-Newtonian model (widespread model 
limitation) with that of inflow rates (major cause of uncertainty). 

Our results demonstrate that TAWSS reacts more sensitively to flow changes 
compared to the linear relationship, with WSS varying by 33-44% for a 25% blood 
flow change. Conversely, changes in OSI with respectto inlet flow variations are more 
complex, exhibiting high variability. The non-Newtonian viscosity model notably 
affects CFD outcomes. TAWSS reduction due to the Modified Cross model matches a 
25% blood flow reduction (-33/-42%). A high sensibility to the rheological properties 
of the fluid is also found in OSI values (-23/-56%). These findings underscore the 
need for patient-specific flow data and non-Newtonian rheology in CFD analysis to 
enhance reliability and reduce variability. 


7.4.2 Relation to Others 


Regarding the influence of inlet blood flow assumptions on WSS metrics, Evju [13] 
also noted a TAWSS decrease of -47% for a 25% reduction in blood flow. Similarly, 
OSI results align with existing literature, exhibiting complex behavior that does not 
correlate with inflow variations [28]. 

Addressing the variability linked to inlet flow rates could be alleviated by incorpo- 
rating patient-specific flows measured via Doppler Ultrasound during clinical prac- 
tice. The complexity deepens concerning the impact of the non-Newtonian model. 
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Steinman's review [6] delves into the non-Newtonian controversy, revealing that 
many CFD analyses overlook non-Newtonian viscosity due to computational costs 
outweighing perceived benefits. However, our work unveils a substantial influence of 
non-Newtonian viscosity on near-wall hemodynamic parameters. Our findings echo 
Mahrous' recent review on non-Newtonian CFD models of intracranial aneurysms, 
demonstrating that the commonly used Newtonian assumption leads to an overesti- 
mation of aneurysm WSS by 17-50% [29]. This effect extends to OSI, as highlighted 
by Oliveira [30], who observed an underestimation of surface-averaged OSI by over 
30%. Thus, we sustain that incorporating a non-Newtonian viscosity model is pivotal 
and should be a consideration in future CFD studies on cerebral aneurysms. 


7.4.3 Limitations 


There are certain aspects of the patient-specific pipeline [6] that have not been 
explored in our discussion, representing limitations of our work and potential avenues 
for future research. For example, we have not taken into account the variability of 
rheological properties among individuals. As per [31], viscosity can vary by up to 
+20% across different subjects. With regard to boundary conditions, our study did 
not explore the effects of different outflow conditions. We applied an inlet blood flow 
and waveform from literature to all models, regardless of parent artery diameters. 
The latter might lead to elevated WSS values in vessels with smaller calibers. 

Further limitations include the lack of an assessment of imaging and lumen 
segmentation, known sources of error [32, 33]. Moreover, our post-processing ex- 
clusively targets the aneurysm sac region, with the selection of the cut zone at the 
aneurysm neck being arbitrary and operator-dependent, thereby introducing vari- 
ability. Numerical considerations pose another limitation, as we did not extensively 
examine convergence in both space and time. However, the adopted spatial and 
temporal discretization guarantees the accuracy of the results. Lastly, our quantita- 
tive analysis focuses on the calculation of TAWSS and OSI, while numerous other 
hemodynamic indicators in the literature could be explored. 
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